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SYNOPSIS 


This  report  presents  the  results  of  a  theoretical  study  of 
inelastic  material  behaviour  produced  by  static  loading  on 
a  viscous  continuum.  A  number  of  issues  are  addressed  in 
an  attempt  to  introduce  inelastic  behaviour  into  piecewise 
numerical  methods.  Methods  based  on  similitude  and  direct 
integration  of  the  displacement  field  are  investigated. 
Potential  field  theory  is  invoked  as  a  means  of  describing 
stress  distributions  which  are  assumed  invariant  with  respect 
to  time  of  straining  the  material.  No  attempt  has  been 
made  to  include  forced  flow,  strain  hardening  or  frictional 
material  properties.  The  analysis  is  restricted  to  slow 
deformations  and  results  in  a  line  integral  for  conversion 
of  strain  rate  to  displacement  as  a  function  of  time  of 
straining.  The  interpretation  of  experimental  data  is 
viewed  from  the  stand-point  of  functional  analysis  as  it 
is  felt  that  this  branch  of  mathematics  has  unexplored 
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CHAPTER  1 


Theoretical  Background  of  Creep  Analysis 


1.1  Introduction 

This  report  covers  the  period  from  August  1976  to 
December  1981  for  research  under  Grant  No.  DA  ERO-76- 
G-o63.  The  object  of  the  research  stipulated  in  the 
grant  documentation  was  to  compliment  earlier  experimen¬ 
tal  study  of  flexible  pavements  by  fundamental  analysis 
of  materials  under  stress.  The  main  objective  is  to 
find  solutions  for  unrestrained  flow  such  as  the  rutting 
of  flexible  payments  and  settlement  of  structures. 

The  original  proposal  was  accepted  6  August  1976  and 
subsequently  was  the  subject  of  one  modification 
designated  DRXSN-E-AO,  of  March  1977. 

This  report  presents  the  results  of  a  theoretical  study 
of  inelastic  material  behaviour  in  a  stress  field 
produced  by  static  loading  of  a  viscous  continuum. 

During  the  course  of  the  project  it  transpired  that  the 
original  objectives  became  somewhat  eclipsed  by  theoretical 
development  of  wider  scope.  Procedures  which  employ 
the  mathematical  background  of  continuum  mechanics  to 
the  solution  of  the  longstanding  engineering  problem  of 
slow  creep  of  materials  are  described.  While  the  theme 
of  application  to  highway  engineering  underscores  the 
study  in  this  report  the  development  of  ideas  is 
generalised  to  encompass  a  variety  of  problems  in  geo¬ 
technical  engineering.  No  attempt  has  been  made  to 
confine  the  study  to  specific  materials,  but  for  the 
sake  of  some  simplification  the  discussion  of  certain 
types  of  inelastic  bahaviour  have  been  omitted  e.g.  flow 
of  granular  materials,  strain  hardening  substances. 
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A  number  of  issues  are  addressed  in  an  attempt  to 
introduce  inelastic  behaviour  into  piecewise  nume  i 
construction  of  solutions  of  some  generality  in  th 
mechanics  of  solids.  The  projec  commenced  ^ith  a 
review  of  existing  theory  and  so'  ition  methou  to  be 
found  in  the  extensive  literature  on  creep  and  plastic 
flow.  The  author's  previous  experience  of  unresolved 
problems  in  geotechnical  engineering  provided  motives 
and  guidance  to  the  ultimate  requirements  of  further 
theoretical  development.  The  text  will  demonstrate 
that  no  definitive  line  is  taken  in  the  preparation  of 
this  report:  rather  the  problem  is  attacked  from  a 
number  of  standpoints  in  a  search  for  common  ground  of 
a  unified  theory;  the  problem  is  treated  at  different 
levels  of  complexity,  but  generally  the  aim  is  simplific¬ 
ation  of  existing  analytical  methods. 

These  methods  comprise  classical  plasticity  which 
exhibits  path  dependent  but  time  independent  behaviour; 
classical  viscoelasticity  and  creep  theory  both  of  which 
have  path  dependent  and  time  dependent  components.  To 
avoid  infringing  on  the  technical  meaning  of  visco¬ 
elastic  and  plastic  behaviour  the  theoretical  developments 
herein  may  be  described  as  that  which  applies  to  material 
in  the  mastic  state  e.g.  asphalt  or  remoulded  clay.  A 
brief  statement  of  the  classical  theories  follows. 

Plasticity  The  task  of  plastic  theory  is  twofold: 
first  to  set  up  relationships  between  stress  and  strain 
which  describes  adequately  the  observed  plastic  deformat¬ 
ion  of  the  material  in  laboratory  tests,  and  second  to 
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develop  techniques  for  using  these  relationships  J.n 
analysis  and  design.  At  least  two  variants  of  the 
theory  are  in  common  usage  -  the  ideal  rigid  plastic  and 
the  elastoplastic  models.  These  models  employ  the 
concepts  of  yield  surface,  associated  flow  rules  and 
plastic  potential  functions.  Throughout  the  theory 
time  enters  as  a  dummy  variable  only,  while  stresses 
and  deformations  are  viewed  as  incremental  rather  than 
total  entities.  There  is  just  one  expression  in 
plasticity  theory  that  conceivably  has  a  bearing  on  the 
creep  problem,  i.e.  the  expression  for  strain  in  terms 
of  plastic  potential  viz. 


de  . 
1 

dt 


=  A 


6f 

6a . 
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where  f  denotes  a  yield  function 
"  a  denotes  stress. 
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The  theory  poses  the  problem  that  for  real  materials  it 
requires  considerable  experimental  work  to  define  the 
yield  function  f(o). 


Viscoelasticity  The  theory  of  linear  viscoelasticity 
is  well  established  and  may  be  summarised  in  the  set  of 
expressions,  Adeyeriand  Krizek  (1969) : 
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Creep  law 
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where  G  and  J  denote  relaxation  moduli  and  creep  com¬ 
pliances  respectively, 

S.  .  denotes  the  deviatoric  stress  tensor 

e .  .  "  the  strain  tensor 

ID 

e.  .  "  the  deviatoric  strain  tensor 

ID 

x  =  (x  ,x  ,x  )  is  a  fixed  Cartesian  co-ordinate  system. 
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The  theory  is  based  on  the  premise  that  stress  depends 
principally  on  the  rate  at  which  the  material  is  deformed 
according  to  the  rheological  relationship 

de 


o  . 
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i 


n- 


dt 


1  =  1 


wnere  E  denotes  Young's  modulus 

n  is  the  viscous  coefficient 


The  theory  is  particularly  applicable  to  problems  of 
forced  flow  where  boundary  stresses  are  clearly  a  function 
of  the  rate  of  deformation  -  as  for  instance  extrusion 
of  a  tube  of  toothpaste.  In  free  flow  problems  the 
theory  is  less  satisfactory  as  in  certain  cases  it 
predicts  deformation  that  tends  to  unrelasitic  values 
as  the  time  scale  is  extended.  This  feature  appears  in 
the  application  of  viscoelasticity  to  rutting  of  flexible 
pavements.  Thrower  (19  77)  . 

Viso-plasticity  This  theory  aims  at  generalisation 
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of  viscous  and  plastic  deformation  on  the  assumption 
that  the  material  behaves  as  an  elastic  solid  exnibiting 
a  zero  rate  of  straining  for  stresses  which  are  below  a 
threshold  or  yield  value.  When  the  threshold  stress 
is  exceeded,  flow  begins  at  a  rate  which  is  a  function 
of  the  excess  stress. 


The  strain  rate  is  related  to  the  yield  function  in  the 
form,  Zienkiewicz  and  Godbole  (1975) 


1.3 


where  p  is  some  viscous  parameter  and  F(a)  represents 
the  yiela  function, 

Hence  the  notation 

<F>  =  0  if  F  <  0 

<  F  >  =  F  if  F>0 


The  resemblance  between  equations  1.1  and  1.3  is  obvious 
and  for  p  =  0  the  visco-plastic  and  ideal  plastic 
formulations  yield  identical  results. 


Theory  of  Creep  The  various  formulations  of  the 
theory  of  creep  are  collected  in  the  works  of  Penny  and 
Marriott  (1971) ,  Rabotnov  (1953) ,  Arutyunyan  (1966)  and 
several  other  authors.  The  treatise  due  to  N.  Kh. 
Arutyunyan  is  of  particular  interest  because  tor  the 
most  part  it  deals  with  analyses  of  unrestrained  creep 
flow  in  civil  engineering  applications.  The  advanced 
mathematical  treatment  undertaken  for  the  creep  of 
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concrete  structures  serves  to  highlight  the  complexity 
of  closed  form  analysis.  However  his  analysis  dis¬ 
closes  very  practical  results  which  are  stated  in  the 
form  of  theorems  thusiy: 

Theorem  1:  If  the  stress  condition  of  a  given  body 

is  produced  by  the  action  of  external  forces,  and  its 
creep  function  for  uniaxial  stress  C(t,x)  is  proport¬ 
ional,  with  constant  coefficient  of  proportionality  kQ 
to  the  creep  function  for  pure  shear  w(t,x)  then  the 
system  of  stresses  in  the  body  considered  coincides 
identically  with  the  system  of  stresses  of  the 
corresponding  instantaneous  elastic  problem. 

Theorem  2:  If  the  stress  condition  in  a  body  either 

is  constant  or  changes  linearly  with  the  x,y,z  co¬ 
ordinates  then  the  stresses  in  the  body  during  creep 
will  coincide  with  the  stresses  of  the  corresponding 
instantaneous  elastic  problem  for  the  same  body  with 
different  coefficients  of  lateral  compression 

v  (t)  and  v  { t , t ) 

1  2 

The  theorems  are  the  outcome  of  the  following 
assumptions : 

1.  The  material  is  regarded  as  a  homogeneous 
isotropic  body; 

2.  The  relationship  between  creep  deformation  and 
stresses  is  linear; 

3.  The  law  of  superposition  applies  to  creep 
deformation . 

4.  The  functions  that  characterize  the  changes  in 

the  coefficient  of  lateral  compression  for  elastic 

deformation  and  creep  deformation  are  identical, 

i.e.  v  (t)  =  v  (t,x)  =  v. 

1  2 
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The  theorems  are  interpreted  to  mean  that  for  un¬ 
restrained  slow  creeping  flow  of  a  material  the  stress 
distribution  is  not  significantly  different  from  that 
of  the  elasticity  analysis  for  steady  state  boundary 
stresses.  Of  course  if  the  boundary  conditions  permit 
a  relaxation  of  stress  anywhere  in  the  domain  then  the 
stress  distribution  will  adjust  itself  to  minimise 
stress  concentration  as  is  the  case  with  forced  flow. 


The  approach  of  Arutyunyan  is  essentially  that  of  the 
hereditary  integral  method  which  originated  with 
Boltzmann  and  was  further  developed  in  the  works  of 
Volterra.  The  fundamental  expressions  of  the  theory 
in  Arutyunyan 's  notation  (see  page  10  of  his  book) 
reads  as  follows: 


exlt) 


ax(t)  [l+v^  (t)]-vi  (t)  (S(t) 
E(t1 


+  (x,y,e)  1.4 

where 

<5  ( t ,  T )  =  +  C<t,x>; 
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V  (-r) 

6  (t,t)  =  — ± -  +  v  (t , t ) C  (t , t)  ; 

1  E(t)  2 

S(t)  =  ox(t)  +  oy(t)  +  ag(t) 

Remaining  relations  obtained  by  cyclic  permutation  of 
x,y,s.  in  above  formulae  all  the  stress  components 
begin  to  act  simultaneously  at  time  x  =  t  .  These 
basic  equations  of  the  hereditary  integral  method 
describe  the  process  of  deformation  in  a  body  by 
taking  into  account  the  changes  of  both  its  volume  and 
its  shape. 

The  difficulty  in  the  construction  of  a  theory  of  creep 
by  this  approach  consists  in  the  choice  of  the  kernels 
in  the  integral  equations  on  the  basis  of  which  solutions 
may  be  obtained  for  the  fundamental  problems  of 
equilibrium  of  a  creep-elastic  medium. 

Arutyunyan  extended  the  linear  theory  to  non-linear 
creep  where  the  constitutive  relationships  do  not 
entail  linearity  between  stress  and  creep  deformation 
at  specific  times.  He  found  that  the  governing 
differential  equation  is  the  generalised  Riccati 
equation  on  the  assumption  that  the  stress-strain- 
time  relationship  is  only  slightly  non  linear.  The 
problem  of  non-linear  creep  was  found  to  reduce  to  the 
solution  of  an  equation  of  the  form: 

^  =  bE  u ( t) [e ' (t)  +  ye(t)l 

dt2  dt  o  l  J 

where  aft)  =  —  and  the  remaining  notation  is 
defined  in  the  original  text  pp.  264-267. 
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1.2  Problem  Statement 

The  foregoing  review  of  theory  serves  to  illustrate 
the  complex  nature  of  any  attempt  at  a  closed  form 
solution  of  material  flow  under  stress.  The  formulation 
of  constitutive  stress-strain  relations  presents  a 
formidable  challenge  and  to  incorporate  these  in  a 
boundary  value  problem  is  even  more  exacting.  Con¬ 
sequently  simplified  numerical  schemes  is  a  desirable 
goal.  The  fact  that  the  geometrical  configuration  is 
not  the  prime  source  of  difficulties  suggests  that  an 
analysis  developed  for  a  particular  geometry  may  be 
readily  extendable  to  other  configurations.  In  this 
work  the  author  will  concentrate  on  so-called  axisymm- 
etric  flow  problems  with  the  conviction  that  the  develop¬ 
ment  could  equally  apply  to  plane  strain  or  plane  stress 
problems  albeit  with  some  modification  to  the  analyses. 
The  axisymmetric  case  is  of  common  occurrence  in 
foundation  engineering  as  for  instance  a  circular 
storage  tank  on  soft  ground,  wheel  loads  on  an  asphalt 
pavement  or  a  circular  cofferdam  on  a  compressible 
stratum. 

The  boundary  value  problem  resembles  the  classical 
Boussinesq  problem  -  that  of  a  concentrated  or  distrib¬ 
uted  load  on  the  surface  of  a  semi- infinite  half-space. 
The  thrust  of  the  analysis  will  be  to  develop  a 
numerical  technique  tor  estimating  the  displacements 
beneath  a  flexible  uniformly  loaded  circular  contact 
area  located  on  the  free  boundary  of  a  semi-infinite 
half  space  of  homogeneous  mastic  material.  The 
analysis  will  delve  into  flow  of  thick  cylinders  subject¬ 
ed  to  radial  forcing  pressures  (internal  and  external) . 
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The  task  of  developing  the  numerical  procedure  consists 
of  the  following  steps: 

i.  postulate  a  mechanism  of  flow  i.e.  establish 
the  geometry  ot  a  plausible  mode  of  flow 
deformation. 

ii.  for  an  incremented  flow  of  this  mechanism 
integrate  the  work  dissipated  in  non-recover- 
able  deformation  over  the  whole  domain  (or 

a  part  thereof) . 

iii.  equate  the  energy  supplied  by  the  forcing 
pressure  to  the  internal  work  and  hence  find 
the  displacement  at  the  source  of  the  disturb¬ 
ing  force  as  a  function  of  elapsed  time  of 
loading . 

iv.  investigate  method  of  substituting  unprocessed 
test  data  for  constitutive  relationships  in 
the  analysis. 

This  scheme  constitutes  the  direct  method  ot  solution 
and  undoubtedly  presents  a  number  of  obstacles,  but  it 
is  in  the  spirit  of  finite  element  analysis.  The 
main  difficulty  is  that  of  specifying  the  constitutive 
relations  for  deformation  in  mutually  perpendicular 
frames  of  reference  i.e.  creep  laws  in  terms  of  stress 
invariants.  On  the  other  hand  an  indirect  method 
based  on  similitude  offers  an  alternative  scheme. 

The  application  of  similitude  and  dimensional  analysis 
is  known  to  produce  practical  solutions  to  otherwise 
intractable  problems,  Glen  Murphy  (1953).  Both 
possibilities  will  be  investigated  in  this  report. 
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CHAPTER  2 

A  Method  of  Calculating  Time-dependent  Deformation 


2.1  Introduction 


A  draw-back  to  the  methods  outlined  in  Chapter  1  is  the 
necessity  to  idealise  the  response  of  real  materials  to 
such  an  extent  that  the  stress-strain  relationship  fits 
into  the  scheme  of  a  particular  theory.  The  constitutive 
equations  tend  to  reflect  the  type  of  proposed  analysis 
rather  than  actual  behaviour  i.e.  rheological  models  for 
viscoelasticity  and  flow  rules  for  plasticity.  Whereas 
the  real  behaviour  may  reflect  a  variety  of  responses, 
the  classical  theories  deal  with  preselected  components. 
This  approach  results  in  a  proliferation  of  material 
parameters,  constants  and  exponents  in  the  effort  to  fit 
experimental  data.  In  what  follows  the  author  proposes 
an  alternative  method  for  the  special  case  of  axisymmetric 
flow  under  the  action  of  sustained  loading.  The  basic 
idea  is  to  relate  the  flow  problem  to  an  experimental 
investigation  of  a  corresponding  model  of  the  stress  state 
using  the  minimum  of  data  processing  of  the  test  results. 
In  principle  the  method  should  prove  equally  valid  for 
both  linear  and  non-linear  behaviour  provided  the  stress 
state  can  be  simulated  in  the  test  procedure.  To  fix 
ideas  we  consider  only  axisymmetric  flow  in  non-forced 
boundary  value  problems;  such  processes  of  flow  as 
metal  forming  with  dies  or  presses  are  excluded  at  this 
stage  in  the  development  of  the  theory. 

Initially  theoretical  aspects  are  explored,  one  of  which 
is  presented  in  this  Chapter.  The  theoretical  exercises 
are  merely  the  forerunners  to  the  main  thrust  of  the  work, 
namely  the  preparation  of  a  scheme  for  numerical  analysis 
by  computer. 
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2.2  Theoretical  Considerations 


The  time-dependent  behaviour  under  consideration  is  a 
slow  process  commonly  known  as  creep.  Because  inertial 
effects  are  excluded  the  laws  of  statics  are  applicable 
in  the  set  of  equipibrium  equations  of  elasticity  theory. 

For  the  purpose  of  this  analysis  it  is  appropriate  to 
consider  the  equilibrium  equations  with  respect  to 
curvilinear  coordinates.  Curvilinear  coordinates  must 
be  considered  as  being  embedded  in  the  material  and  are 
defined  in  terms  of  a  function  which  is  assigned  some 
value  at  each  point  throughout  the  material.  The  direction 
and  curvature  of  the  coordinates  therefore  changes  from 
point  to  point.  Let  the  first  family  of  curves  be 
defined  by 

f(x,y)  =  a 

and  the  second  family  be  defined  by 
g(x,y)  «  3 

where  x  and  y  are  rectangular  coordinate  components. 

For  the  case  of  a  two-dimensional  stress  state  referred 
to  curvilinear  coordinates  the  equilibrium  equations 
are  given  by  the  expressions.  Ford  and  Alexander  (1976) s 
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Where  p„  and  pQ  denote  radii  of  curvature,  and  o  and  t 
a  p 

are  direct  and  shearing  stresses,  respectively. 
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If  the  coordinate  axes  are  chosen  to  coincide  with  the 
directions  of  principal  velocity  components,  then  for 
materials  that  exhibit  sliding  on  planes  inclined  at  45 
degrees  to  the  principal  planes,  the  equations  of 
equilibrium  in  two  dimensions  along  lines  of  maximum 
shear  stress  can  be  written  in  the  form,  Kuske  and 
Robertson  (1974) 
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where  o  ,o  are 
1  2 

p  ,p  denote  the 
1  2 

maximum  shear 


the  principal  stresses  at  the  point  and 
radii  of  curvature  of  the  lines  of 


FIG.  2.1: 


Element  Bounded  by  Lines  of  Maximum  Shear 
Stress 
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Equation  2.2  provides  the  equilibrium  conditions  for 
the  state  of  stress  shown  in  Fig.  2.1.  This  state 
consists  of  the  distributions  of  the  mean  normal  stress 
and  shears  along  the  trajectories  of  the  curvilinear 
coordinate  system.  The  stresses  at  any  point  within  a 
cell  reduce  to  the  mean  normal  stress  and  the  maximum 
shear  stresses.  If  we  confine  the  analysis  herein  to 
materials  that  flow  along  planes  of  maximum  shear  stress 
the  net  is  geometrically  an  orthogonal  mesh.  Materials 
that  possess  a  Coulomb  frictional  component  are  thus 
excluded.  In  the  notation  of  the  slip-line  theory  of 
plasticity  the  geometry  is  termed  an  ot , 3  net.  The 
a, 3  lines  can  be  assigned  curvilinear  coordinate  values 
according  to  scale  of  plot  but  one  of  the  lines  -  the 
3  line  -  has  the  physical  significance  that  it  is  the 
contour  of  constant  normal  stress.  The  3  -  lines  are 
readily  determined  as  the  contours  of  mean  stress  for 
boundary  value  problems  by  solution  of  Laplace's 
equation  viz. 
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2.3 


which  follows  from  the  invariance  of  the  Laplacian 
under  transformation  of  coordinate  axes.  The  equation 

2.3  provides  the  basis  for  the  photoelastic  technique 
because  it  implies  that  stress  distribution  is  independent 
of  the  elastic  constants  of  a  homogeneous  medium.  The 
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conjugate  function  gives  the  velocity  characteristics 
a  -  lines*.  It  would  be  a  simple  matter  to  plot  the 
locations  of  the  -  lines  if  the  shear  stress 
component  could  be  used  as  an  argument  of  Laplace's 
equation,  but  since  this  is  not  the  case  harmonic 
analysis  only  partially  solves  the  problem  of  plotting 
the  slip-line  field.  The  method  of  stream  functions 
has  been  used  to  find  the  geometry  of  the  a  -  lines 
separately  from  the  0  -  lines  by  many  investigators 
in  the  field  of  fluid  mechanics.  In  Chapter  5  of 
this  report,  a  method  of  locating  the  a  -  lines 
is  presented. 


*  The  term  velocity  characteristic  is  borrowed  from 
plasticity  theory  notwithstanding  the  fact  that  the 
governing  equation  2.3  is  elliptic. 
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2.3  A  Theory  of  Flow  based,  on  Stress  Gradient 


Consider  two-dimensional  stress  state  with  reference 
to  curvilinear  coordinates  as  depicted  in  Fig.  2.2. 
Inspection  of  the  Equations  2.1  &  2.2  reveals  that  the 
source  of  viscous  flow  is  the  activation  produced  by 
the  first  term,  namely,  the  gradient  of  the  mean  stress 
normal  to  the  direction  of  flow,  i.e.  the  mean  normal 
stress  variation  alona  a  -  lines  of  Fig.  2.2. 
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FIG.  2.2:  Stress  Distribution  in  Cell  of  the  a-&  Net 
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Equilibrium  is  preserved  by  the  shear  stress  distribution 
on  intersecting  trajectories  as  shown  by  the  remaining 
terms  in  the  equations.  In  particular  for  large  radii 
of  curvature  of  the  trajectories  the  third  term  of 
Equation  2.2  makes  only  a  small  contribution  to  equilibrium. 
Now  in  a  viscous  continuum  the  shearing  stresses  are  not 
able  to  completely  restrain  permanent  deformation  once 
the  material  has  absorbed  the  elastic  strain  energy  of 
the  stress  state.  Thus  the  phenomenon  of  creep  ensues 
or  for  certain  boundary  conditions  a  relaxation  of  stress 
occurs. 

The  stress  gradient  along  a  flow  path  bounded  by  two 
a  -  lines  is  induced  by  the  boundary  traction  because 
the  source  is  the  external  pressure  on  the  zone  bounded 
by  the  same  a  -  lines.  The  task  is  to  relate  the  creep 
in  a  model  test  to  the  stress  gradient  in  the  test  and 
subsequently  to  use  the  correlation  for  the  purpose  of 
predicting  flow  in  a  prototype  of  the  boundary  value 
problem.  To  accomplish  the  transition  from  model  to 
prototype  it  is  essential  that  a  unique  expression  of 
the  stress  gradient  be  found.  If  the  material  behaviour 
is  linear  the  so-called  unique  expression  need  only  be 
determined  for  a  single  value  of  the  stress  state,  but 
for  non-linear  material  a  relationship  over  a  range  of 
values  of  stress  presumably  is  required.  Such  a 
unique  relationship  will  advance  the  theory  to  the 
extent  that  it  provides  the  link  between  driving  forces 
in  model  and  prototype. 


(17) 


In  the  search  for  a  unique  correlating  function  the 
author  was  attracted  to  the  properties  of  the  Laplace 
transform.  It  is  well  known  that  the  transform  confers 
the  capability  of  distinguishing  between  functions  which 
produce  the  same  value  for  their  integrals  over  some 
one  interval.  By  definition  we  have  the  Laplace 
transform  with  independent  parameter  p. 


L{f(x)}  =  j“szof(x)e~FXdx  2.4 

Clearly  it  is  not  appropriate  to  assign  values  to  p 
if  we  are  searching  for  a  unique  expression  in  the 
hope  that  a  function  such  as  the  distribution  of  stress 
gradient  can  be  characterised  by  a  single  numerical 
value.  These  considerations  have  led  to  the  intro¬ 
duction  of  a  special  function  which  bears  resemblance 
to  Laplace  transformation  formulae;  it  has  the  same 
property  -  that  of  attenuating  the  function  towards  the 
end  of  its  domain.  The  proposed  function  has  the  form 
_  x-a 

e  b  x  which  replaces  e  px  of  the  Laplace  transform. 

By  analogy  we  can  write  for  a  flow  path  the  unique 
expression  for  distribution  of  f(x)  as  follows: 


x-a 

G{f(x)}  =  ~Xdx 

[a,b] 


2.5 


where  a  and  b  are  the  extremal  points  of  path  p. 
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The  correspondence  between  prototype  and  model  can  now 
be  postulated  by  assuming  that  for  equal  values  of  G 
and  the  same  area  subjected  to  surface  traction  the 
power  expended  in  producing  permanent  deformation  is 
identical;  and  where  the  G  function  is  different  for 
prototype  and  model  the  boundary  velocities  for  the 
prototype  follow  from  similitude  according  to  the 
relationship: 
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where  a  denotes  the  driving  stress  on  the  boundary 
v  denotes  velocity  at  the  source 
and  subscripts  p  and  m  mean  prototype  and  model, 
respectively. 


To  illustrate  the  physical  meaning  of  Equation  2.6  we 
take  for  our  model  the  radial  flow  in  the  wall  of  a 
thick  cylinder  under  internal  pressure  and  for  the 
prototype  the  axially  symmetric  deformation  of  a  semi¬ 
infinite  viscous  continum  subjected  to  uniform  pressure 
on  a  circular  contact  area  as  shown  in  Fig.  2.3-  To 
achieve  a  correspondence  we  must  compare  the  G-values 
on  paths  that  have  corresponding  stress  states. 

Firstly  consider  the  case  of  linear  material  behaviour. 
An  expansion  test  on  a  thick  cylinder  will  provide  the 
information  to  solve  Equation  2.6  because  the  only  un¬ 
known  quantity  is  the  velocity  (u^)  ;  the  velocity  will 
in  general  be  a  function  of  time  and  hence  displacements 
will  need  to  be  determined  by  integration  over  an 
interval  of  time  or  alternatively  as  a  set  of  Riemann 
sums.  Second  for  non-linear  materials  the  only 


(20) 


departure  from  the  procedures  outlined  above  is  that 

the  relationship  between  the  velocity  and  the  transform 

G  be  determined  by  test  on  a  suitable  model  of  the 

stress  state  over  a  range  of  the  stresses.  The  whole 

process  is  considerably  simplified  if  the  assumption 

of  incompressibility  of  the  material  is  tenable.  For 

many  materials  susceptible  to  creep  this  is  approximately 

true  once  the  elastic  deformations  have  developed. 

The  term  in  t^ie  closed  interval  (a,b)  produces  the 

e 

decay  expression  with  numerical  values  as  follows:- 
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The  scheme  can  now  be  applied  to  our  problem  as 
depicted  in  Fig.  2.4.  The  velocity  field  is  taken 
along  a  line  midway  between  two  a  -  lines  which  gives 
the  path  of  integration  of  the  function  G{ai(a)} 

p 

where  o^(a)  is  the  gradient  of  the  normal  stress  in 
the  a-direction.  The  forcing  function  is  simply  the 
boundary  pressure  multiplied  by  the  area  of  the  flow 
'tube'  at  the  boundary. 


FIG.  2.4:  Axi symmetric  Flow  -  an  a,  3  net  with 

Integration  Path 


The  numerical  value  of  G  for  the  flow  tube  ABCD  is 
the  path  integral  along  p(ABCD)(viz. 
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where  a  is  the  distance  along  the  flow  path  measured 
from  the  source  at  aQ  and  afc  is  the  terminal  point  on 
path  p. 

The  flow  tube  cannot  be  taken  in  isolation  because  it 
is  bounded  by  contiguous  zones  of  the  a, 0  net  which 
are  also  in  a  state  of  compatible  flow  i.e.  there  are 
no  discontinuities  in  the  velocity  across  a  -  lines. 


FIG.  2.5: 


Radial  Flow  in  Hollow  Cylinder  Test 


2.4  Determination  of  the  Displacement  Field 

If  a  material  is  un compressible  the  volume  occupied  is 
constant;  the  quantity  of  material  which  enters  a  sub- 
region  is  equal  to  that  which  leaves  it  according  to 
the  continuity  equation  of  fluid  mechanics 
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The  movement  of  the  material  can  be  described  in 
either  Eularian  or  Lagrangian  coordinates;  in  any 
case  the  relationship  between  the  two  systems  is  well 
established,  Hodge  pp.  143-144. 


Due  to  the  comparatively  small  scale  of  movements  in 
a  highly  viscous  medium  the  Lagrangian  approach  is 
deemed  adequate  for  tracking  the  progress  of  flow  from 
one  zone  to  another;  the  change  in  stress  gradiant  is 
sufficiently  small  to  warrant  the  adoption  of  a  fixed 
reference  frame.  Thus  we  will  focus  on  the  relationship 
that  expresses  the  displacements  caused  by  flow  of 
material  into  and  out  of  a  cell  of  a  fixed  a,fc(  net. 

The  displacements  are  to  be  determined  along  the  median 
path  which  is  midway  between  adjacent  a  -  lines,  i.e. 
the  path  p  shown  Fig.  2.6. 

The  method  proposed  is  to  use  the  properties  of  the 
Jacobian  functional  derivative  for  comparison  of  areas 
in  the  two-dimensional  net  and  the  relationship  between 
curvilinear  and  rectangular  coordinates  for  the  third 
dimension  of  the  flow  tube.  The  constant  volume 
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condition  can  be  expressed  in  terras  of  the  distribution 
of  the  Jacobiam  along  the  flow  path  and  the  variation 
in  area  of  the  flow  tube  normal  to  direction  of  flow. 
With  reference  to  Fig.  2.6  the  relationship  in  terms 
of  displacements  is  derived  as  follows: 

Let  J  denote  the  Jacobian  defined  by* 


6a 

6a 

6a; 

6(o,6) 

63 

6a: 

66 

6  ix,y) 

The  numerical  value  of  the  Jacobian  varies  along  the 
median  and  it  is  a  regular  function  (non  vanishing)  as 
demonstrated  in  the  next  Chapter.  Thus  the  displace¬ 
ment  6^  normal  to  the  6  -  lines  is  given  by  the 


FIG:  2.6  Sketch  of  Flow  Tube 


*  The  topic  of  image  mapping  via  the  Jacobian  transform¬ 
ation  is  discussed  in  Chapter  3  of  this  report.  The 
family  of  curves  f(x,y)  =  n  and  g(x,y)  =  6  are  consider 
ed  determinable  entities  at  this  stage. 
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expression 
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The  displacement  of  the  centroid  of  a  cell  6s  is  given 
by  the  linear  approximation: 


6s 


-  +  6fW 


2.11 


where  the  notation  is  that  shown  in  Fig.  2.6.  On  the 

boundary  the  displacement  in  a  given  interval  of  time 

due  to  surface  traction  reduces  6  to  60  if  we  take 

s  j  Pl 

the  first  cell  as  a  line  element 


i.e.  n  =  1  6b  =  6s 
i  i 

where  6s  is  deduced  from  the  velocity  as  given  by 
Equation  2.6.  The  Equations  2.10  and  2.11  enable 
tracking  the  values  of  displacements  from  the  source 
of  boundary  perturbation  to  a  terminal  point  on  the 
boundary  i.e.  the  exit  of  the  material  from  the  net. 


To  calculate  the  displacement  it  is  necessary  to  have 
an  expression  for  the  cross  sectional  area  of  the  flow 
tube  at  stations  along  the  flow  path.  The  relation¬ 
ships  of  the  theory  of  curvilinear  coordinates  provide 
a  general  expression  for  the  area,  Marsden  and  Tromba 
(1975).  pp  3i4. 


viz . 
where 
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Hz)  denotes  the  area  at  any  section 
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2.5  Synthesis  of  Theoretical  Formulations 
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The  foregoing  analysis  involved  the  hypothesis  that 
the  rate  of  flow  can  be  related  to  the  gradient  of 
mean  normal  stress  a  priori.  The  mean  normal  stresses 
are  considered  spatially  invariant  with  time  and 
those  values  determined  by  linear  elastic  theory  for 
constant  boundary  conditions.  The  notions  of 
similitude  and  cross-correlation  were  applied  to  yield 
a  dimensionally  consistant  relationship,  i.e.  equation 
2.6.  In  order  to  establish  a  connection  between  the 
author’s  approach  and  orthodox  analysis  of  flow  problems, 
the  governing  differential  equations  for  viscous  flow 
is  invoked.  These  classical  equations  collectively 
known  as  the  Navier-Stokes  equations  are  basic  to 
analysis  to  both  Newtonian  and  non-ideal  fluids. 


The  Navier-Stokes  equation  for  two-dimensional  steady 
flow  of  an  incompressible  fluid  with  constant  coefficient 
of  viscosity  n  and  density  p  reads: 
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2.13(a) 


Where  u  and  v  are  the  velocity  components  in  the 
x,y  coordinate  directions  and  p  is  the  mean  fluid 
pressure  at  a  point  (not  a  hydrostatic  pressure) 

Chang  Lu  (1973).  Equation  2.13(a)  together  with  the 
continuity  equation  completes  the  classical  formulation. 
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In  the  two-dimensional  region  the  continuity  equation 
is  given  by  the  expression: 


6  u  6v_ 

i>x  &y 


=  0 


Neglicting  acceleration  but  including  inertial  terms 
we  have  the  expressions 
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2.13(b) 

Under  a  coordinate  transformation  to  another  orthogonal 
frame  of  reference  the  property  of  invariance  enables 
writing  the  equation  2.13(b)  for  the  a, 3  coordinates 
in  the  form: 
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2.13(c) 


where  now  u  and  v  are  referred  to  the  orthogonal  spatial 
coordinates  of  the  a, 6  net,  and  the  term  is  known 
from  a  separate  analysis  of  the  stress  field. 


The  first  equation  of  the  set  2.13(c)  only  need  be 
considered  because  the  mean  pressure  gradient  is  zero 
everywhere  within  the  solution  domain  for  the  right 
hand  side  of  the  second  equation,  i.e. 
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This  nonlinear  and  non-homogeneous  ordinary  different¬ 
ial  equation  can  be  solved  for  u  if  a  particular 
solution  can  be  deduced,  Courant  and  John  (1974) 
Hidlebrand  (1957) .  The  general  solution  is  formed 
from  the  known  function  by  exponentiation  and  the 
ordinary  process  of  integration  with  one  arbitrary 
constant  as  given  in  Courant  1974.  It  appears  that 
such  solutions  are  mainly  relevant  in  fluid  mechanics 
or  aravitational  flow  problems  i.e.  problems  where 
the  driving  forces  are  derived  from  the  gravitational 
influence  on  the  mass  density.  However,  it  is 
significant  to  note  that  the  displacements  will  have 
similar  distributions  whether  equations  2.11  or  2.13(e) 
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are  used.  Finally  the  equation  2.13(e)  has  as  its 
domain  of  independent  variable  a  the  path  of  integration 
inscribed  on  the  physical  flow  line  which  is  the  same 
median  path  as  that  proposed  for  evaluating  the 
correlation  function  Gfo1}  of  equation  2.7. 
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CHAPTER  3 

Mapping  Functions  for  Orthogonal  Nets 
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3  .  i  /•’  i.-)W  nets:  Approx  i  ma  tion  b  y  Po  lynomial  s 

In  this  chapter  a  method  for  the  functional  representation 
of  the  geometry  of  orthogonal  trajectories  is  presented. 
Families  of  orthogonal  trajectories  are  encountered  in 
hydraulic  flow  nets,  electrostatics  and  stress  fields, 
or  indeed  any  physical  problem  described  by  Laplace's 
Equation.  Flow  nets  are  conventionally  derived  from  the 
theory  of  complex  variables  in  an  exact  formulation, 
approximated  by  the  finite  difference  and  finite 
elements  techniques,  or  simply  sketched  freehand  to  match 
a  particular  set  of  boundary  conditions.  The  generating 
functions  deduced  from  complex  variable  theory  can  be 
quite  complicated,  while  no  functional  representation 
emerges  from  the  alternative  methods.  The  problem  as 
posed  herein  is  to  fit  simple  expressions  which  adequately 
describe  the  geometry  in  functional  relationships. 

The  method  proposed  is  that  of  trial  fitting  of  poly¬ 
nomials  where  the  selected  polynomials  consist  of  a  set  of 
harmonic  mapping  functions.  These  functions  are 
generated  from  the  definition  of  a  complex  function  viz. 

f(x,y)  =  (x  +  iy)n  :  n  =  1,2,3  ...;  i  =  /-T 

Thus  by  taking  values  of  n  to  order  4  we  obtain  harmonic 
functions  which  are  deemed  to  approximate  the  geometry  of 
an  orthogonal  net.  Letting  4>  and  tp  represent  the  co¬ 
ordinates  of  a  mapping  function  Ln  a  Cartesian  frame  of 
reference  then  the  geometry  is  described  by  the 
relations : 
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41  =  A(2x+y)  +  B(x-y1 2)  +  C(x3-3xy2)  +  D  (xH+y 4 -6xzy2 )  +  $ 

^  =  A(2y-x)  +  B(2xy)  +  C(3x2y-y3)  +  D(4x3y-4y3x)  +  4*0 

(3.1) 

where  A,  B,  C,  D  are  scaling  factors  (real  numbers) . 

The  functions  <j>  (  x,y)  ,  (x,y)  map  a  rectangular  mesh  into 

a  distorted  shape  where  the  lines  x  =  c^,  y  =  c2  in  the 
(x,y)  plane  appear  as  orthogonal  trajectories  in  the 
(4>,4>)  space  irrespective  of  the  values  of  the  sealing 
factors  or  the  constants  c^,  C2*  Now  suppose  that  a 
plot  of  <f>  versus  ip  in  Cartesian  coordinates  is  available 
from  some  source  problem  in  electrostatics  or 
ground  water  seepage.  The  task  is  to 
firstly  find  the  origin  of  the  generating  functions 
<t>  (x,y)  /  ^  (x,y)  if  in  fact  the  plot  stems  from  a  single 
origin,  and  second  to  extract  the  expressions  that  express 
the  geometrical  pattern  in  algebraic  form.  Assuming 
then  that  the  net  can  be  traced  to  a  single  origin,  i.e. 
not  bipolar  or  multipolar  coordinate  systems,  the  proposed 
analysis  proceeds  as  follows: 

1.  A  cell  in  (4>,i|j)  space  bounded  by  four  intersecting 
trajectories  is  selected  (remote  from  the  confluence 
of  lines  that  indicate  a  possible  location  of  this 
origin  of  the  coordinate  system) . 

2.  Simultaneous  equations  are  written  for  the  selected 

order  of  the  harmonic  functions  by  taking  the  nodes 

in  turn  with  one  of  the  four  corner  nodes  as  a 
temporary  origin. 

3.  The  set  of  simultaneous  equations  is  solved  for  the 


unknown  values  of  the  constants  as  given  in  Eqn.  (3.1) . 

4.  A  search  for  those  values  of  the  constants  which  are 
consistent  with  the  constraint  that  the  coordinate 
intervals  of  (x,y)  are  restricted  to  integers  is 
implemented. 

5.  The  global  origin  with  respect  to  the  temporary  origin 
is  located  by  referring  back  to  the  results  of  the 
search  of  step  4  above . 

Even  in  the  case  of  a  bipolar  net  the  method  can  yield  the 
functional  relationships  (at  least  locally)  and  it  also 
may  be  applied  to  freehand  plots  in  a  piecewise  scan  of 
the  net  cells.  Even  freehand  plots  lead  to  well  con¬ 
ditioned  equations  i.e.  small  changes  in  <p,  ip  values 
are  not  grossly  reflected  in  the  values  of  the  scaling 
factors . 

The  algebra  can  be  considerably  reduced  by  eliminating 
biquadratic  terms,  i.e.  letting  D  =  0;  this  is  a 
simplification  that  is  well  justified  because  of  the 
ability  of  cubics  to  produce  the  optimum  fit  to  an 
arbitrary  continuous  function  (c.f.  the  theory  of  splines 
in  succeeding  chapter) .  Consider  the  plot  of  two 
functions  <p  and  in  Cartesian  coordinates  as  shown  in 
Fig.  3.1.  If  we  isolate  any  cell  and  select  a 
temporary  origin  e.g.  cell  and  the  node  marked  1  of  inset 
in  Fig.  3.1  we  can  write  the  coordinates  of  the  nodes 
as 

•Node  1  (xq,  y0> 

2  *o  +  1] 

3  (x0  +  +  X) 

4  (xo  +  X'  Yo] 
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where  (xQ,yo)  are  the  specific  Cartesian  coordinates  of 
the  lower  left  hand  node  of  the  cell. 

Substituting  these  ordered  pairs  in  Equation  (3.1)  for 
the  abbreviated  cubic  version  we  obtain  the  following: 

Node  1 

*1  “  A(2xo+yo)  +  B(x2o‘y2o)  +  C(x3o’3xoy2o)  +  *o 
=  A<2yo~xo)  +  B(2xoyo}  +  C(3x2oyo’y3o)  + 

Node  2 

*2  =  A(2xo+yo+i)+B(xo-yo_2yo"1)+C(xo"3xoy2- 

6xoyo-3xo}  +  *o 

*2  =  A(2yo~xo+2)+B(2xo+2xoyo)+C(xo+3xoyo‘ 

yo“3yo'3yo_1)  +  ^o  (3  2) 

and  similarly  for  nodes  3  and  4.  Knowing  the  values  of 
<})  and  at  the  nodes  we  can  take  any  three  of  the  above 
equations  and  solve  for  the  constants  A,  B,  C  over  a 
range  of  integer  values  of  X0»Y0*  transpires  that 

if  there  ip  a  unique  set  of  values  of  the  constants  that 
satisfy  all  combinations  of  the  equations  for  the  nodal 
values  throughout  the  net  then  a  global  origin  exists 
for  the  particular  net.  The  temporary  origin  at  the 
unique  values  of  x0fY0  gives  the  clue  to  the  location  of 
global  origin.  The  offsets  are  then  <pQ  and  i|>  in  Equation  3.2. 

Before  a  numerical  example  is  presented  it  is  worthwhile 
to  take  a  look  at  closed  mathematical  solutions  of  the 
problem  of  inverting  images.  With  the  assistance  of  a 
mathematician,  the  author  arrived  at  a  specimen  solution 
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of  a  comparatively  simple  mapping  problem.  The  exercise 
was  motivated  by  remarks  in  a  textbook  to  the  effect  that 
it  should  always  be  possible  to  invert  a  mapping  of 
continuous  functions  (according  to  Riemann’s  Theorem) 
given  the  image Courant  and  John  (1974). 

The  mapping  selected  is  given  by  the  harmonic  functions 

$  =  x2  +  4x  -  y2  +  2y 

=  2xy  +  4y  -  2x 

To  find  the  inverse  we  must  be  able  to  find 
X  =  f1(4>,4')/  y  =  f2(<)>/4') 

Now  4>  =  (x+2)2  -  (y-1)  2-3 

or  <J>+3  =  X2  -  Y2  where  X  *  (x+2),  Y  =  (y-1) 

ip  =  2y  (x+2 )  -  2  (x+2) +4 
=  XY .  whence  Y  =  ^ 

X2(<p+3)  =  x“-(i~)2 

x*-(*+3)x*  •  (^)2  =  0 


and 


or 


iz± 

2 


or 


It  should  be  noted  that  while  Riemann's  mapping 
theorem  demonstrates  the  existance  of  a  mapping 
function,  it  does  not  actually  produce  the  function. 
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X2 


X 


y 


<j>+3  +  /(<J>+3)  2+  C4#-4 )  2 

=  (x+2) 2 

'  -  2  =  f1(*#*) 

+  1  "  £2(*'*>  3'3<b)  °-E‘D 


♦+3  +  /($+3  )z+ 


2 


It  is  instructive  to  refer  to  Courant  and  John's  text 
(or  similar  treatise)  to  see  just  how  difficult  the 
problem  can  become  for  the  general  case  of  image  mapping, 
and  indeed  a  numerical  procedure  is  proposed  as  the  most 
appropriate  method.*  Even  in  the  case  of  the  simple 
mapping  investigated  above  the  quadratic  form  of  the 
inverse  image  precludes  finding  a  unique  mapping  in 
the  ( x , y )  plane.  Indeed  the  effort  required  to  extend 
the  closed  form  approach  to  harmonic  functions  involving 
cubics  is  considerable  as  evidenced  in  attempts  to  apply 
Descartes  method  for  solving  cubic  polynomial  equations. 
Fortunately  the  reverse  image  is  ameniable  to  numerical 
analysis.  The  numerical  technique  is  easily  programmed 
and  although  it  would  be  impracticable  to  present  the 
computer  output  in  full  abridged  results  serve  the  purpose 
of  highlighting  the  main  features  of  the  numerical  scheme. 
In  what  follows  two  typical  cells  are  selected  from  a 
flow  net,  the  net  is  shown  in  Fig.  3.2,  and  the  values  of 
<p  and  tjj  are  read  from  the  geometry  of  the  net  referred 
to  Cartesian  coordinates. 


*  The  numerical  method  merely  enables  mapping  in  the 
neighbourhood  of  an  isolated  point  and  hence  it  is  of 
limited  usefulness. 
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3.2  Example  of  Decoded  Flow  Net 

Suppose  that  part  of  a  net  which  does  not  include  an 
origin  (source  or  sink)  is  given  in  rectangular  co¬ 
ordinates  (x,y) .  The  problem  posed  is  one  of  extract¬ 
ing  the  analytic  function  that  represents  the  entire 
region  of  the  net  i.e.  to  find  the  function 

5  =  <}>  (x,  y)  +  ii|*(x,y) 

where  4>  and  ip  are  assumed  continuously  differentiable 
functions  (excepting  the  neighbourhood  of  the  origin) . 

A  successful  outcome  of  the  exercise  will  enable 
plotting  the  net  over  any  region  of  interest  and  it 
will  also  provide  the  Jacobian  in  algebraic  form. 

A  condition  for  computational  success  is  that  the  mesh 
is  uniformly  divergent  from  its  origin.  The  analysis 
of  this  problem  by  the  method  of  polynomial  harmonic 
functions  is  demonstrated  with  the  aid  of  a  computer 
generated  mesh  as  shown  in  Fig.  3.2.  This  net  was 
plotted  by  assigning  values  to  the  constants  a,B,C  in 
Equation  3.2.  Although  the  answer  is  known  before¬ 
hand  the  main  features  of  the  approach  is  illustrated 
by  selecting  a  precisely  drawn  net. 

The  known  values  of  <p  and  tj;  provide  eight  algebraic 
equations  per  cell.  These  equations  will  in  general 
be  non-linear  in  terms  of  (x,y)  coordinates.  However 
as  the  (x,y)  values  can  be  treated  as  parameters, 
each  taking  on  integer  values  only,  the  Equations  3.2 
can  be  programmed  within  a  loop  over  a  range  of  the 
integers.  Such  a  search  will  yield  the  values  of 
(x,y)'  relative  th  a  common  origin  for  which  the  constants 
A,  B,  C  and  D  are  invariant  with  respect  to  the  position 
of  the  cell. 


1 1 1  ir '~nr>  ~  i  ' - 
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Cubic  and  quartic  interpolation  schemes  are  handled  with 
equal  ease,  so  the  search  for  a  unique  set  of  constants 
applicable  to  each  of  the  cells  can  commence  with  a  cubic 
polynomial  approximation.  The  mapping  is  illustrated  in 
Fig .  3.3. 

Thus  we  can  write  sets  of  simultaneous  equations  for  the 
nodes  1  to  4  as  follows: 


3 


FIG.  3.3  Typical  cell  of  Orthogonal  net 


SET  I 

4>i  =  A(2x+y)  +  B(x2_y2)  +  C(x3-3xy2) 

4>i  =  A (2y-x)  +  B (2xy)  +  C(3xzy-y3) 

<K  =  A(2x+y+l)  +  B(xz-yz~2y-l)  +  C  (3x2y+3x2-y 3-3y-l) 

SET  II 

4> 2  =  A(2x+y+2)  +  B  (xz+2x+l-yz )  +  C  (x3+3x2  +  3x+l-3xy2-3y2 ) 
\p2  =  A(2y-x-l)  +  B(2xy+2y)  +  C  (3xzy+6xy+3y-y 3 ) 

^4  =  A(2y+2-x)  +  B(2xy+2x)  +  C (3x2y+3x2-y3-3ya-3y-l) 

SET  III 

4>4  =  A(2x+y+l)  +  B  (x2-y  z-2y-l)  +  C  (x3-3xy2-6xy-3x) 

4) 4  =  A(2y+2-x)  +  B(2xy+2x)  +  C  (3x2  +  3x2-y s-3y2-3y-l) 

4)2  *  A  (2x+y+2)  +  B(x2+2x+l-y2)  +  C  (x3+3x2+3x+l-3xy2-3y2 ) 
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SET  IV 

4> 2  =  A(2x+y+2)  +  B  (x2+2x+l-y2 )  +  C  (x3+3x2+3x+l-3x+l-3xy2-3y2 ) 
<p  i  =  A(2x+y)  +  B(x2-y2)  +  C(x3-3xy2) 

04  =  A(2x+y+l)  +  B (x2-y2-2y-l)  +  C (x3-3xy 2-6xy-3x) 


SET  V 

02  =  A(2y-x-l)  +  B{2xy+2y)  +  C ( 3x2y+6xy+3y-y 3 ) 

0i  =  A(2y-x)  +  B(2xy)  +  C(3x2y-y3) 

04  =  A(2y+2-x)  +  B(2xy+2x)  +  C (3x2y+3x2-y3-3y2-3y-l) 


The  equations  above  are  solved  for  the  constants  A,  B  and  C 
by  inserting  values  of  0  and  0  at  the  nodes  and  a  range  of 
integer  values  for  x  and  y;  the  process  is  continued  until  a 
particular  pair  of  x,y  values  yield  the  same  values  for  the 
constants  irrespective  of  the  ordering  of  the  set  of  equations. 
The  consistent  values  are  identified  in  Table  111-1 by  arrows  in 
right  hand  margin.  Because  of  re-ordering  of  node  numbers  in 
programs  run  on  the  Hewlett  Packard  and  Digital  computers  the 
diagrams  are  at  some  variance  with  the  tabulated  results. 
Nevertheless,  the  results  as  listed,  the  program  in  BASIC, 
demonstrates  the  key  features  of  the  analysis. 
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Cell  i  SET  I  SET  II  SET  III 


The  results,  determined  for  cells  i  and  j,  are  given  in 
Table  III-l  where  the  constants  for  a  partial  list  of  the 
independent  parameters  are  tabulated. 


The  position  Of  the  global  origin  is  at  the  same  location 
for  the  two  cells  chosen  here: 

Cell  i  :  origin  :  x  =  5  and  y  =  5 
Cell  2  :  origin  :  x  =  7  and  y  =  6 

Constants  A  =  1.00,  B  =  -1.41,  and  C  =  -1.41,  D  =  0. 

The  results  can  be  verified  for  all  cells  outside  the 
immediate  vicinity  of  the  origin  i.e.  where  the  change  is 
curvature  is  moderate.  The  equation  of  the  entire  net  are 
now  established  in  the  relationships 

=  =2x+y-l . 41 (x2-y2 ) -1 . 41 (x3-3xy2 ) 
ip  =  2y-x-2 . 82  (xy) -1. 41  (  3x2y-y 3) 

and  the  analytical  function  ?  is  completely  defined.  The 
flow  net  can  be  extended  indefinitely  by  using  these 
functions . 

However  in  stress  analysis  the  relationship  is  normally 
encountered  with  form  ip  =  f  ( 4> )  .  Thus  equation  3.3  would 
be  reduced  to: 

^  =  2X  (X2-<t>-3)  ^  +  4  (X  =  x+2 )  3.3.C 

On  a  Cartesian  plot  the  expression  i|>  =  f  (<}>),  for  a  specific 
domain  of  4>(x,y),  is  a  family  of  curves  because  i|i  is  a 
multivalued  function  if  we  vary  x  or  y.  The  curvilinear 
plot  is  more  readily  described  in  curvilinear  coordinates 
wherein  the  points  set  (x,y)  are  regarded  as  curvilinear 
coordinates.  From  this  viewpoint  the  differential  relation¬ 
ship  for  areas  of  cells  is  given  by  the  expression 
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d<j>d4  =  |j|dx  dy 

where  |j|  is  the  determinant  of  the  Jacobian  matrix. 


The  value  of  the  Jacobian  is  given  by  the  expression 

|J|  «  t*  .  -  t*  .  t*  3.4 


6  4>  64*  _  <5  <t>  64 

<J>x  6y  6y  6x 


Hence  for  cubic  interpolation 

|j|  =  5A2  +8ABx+12ACx2-12 ACy 2  +4B2x2  +4B2y 2-4ABy+ 

12BCx  3+12BCxy2+9C2xH+18C2x2y2+9C2y‘’-12ACxy 


3.5 


For  interpolation  up  to  and  including  quartic  terms  (n=4) 
the  Jacobean  is  given  by  the  following  terms  inserted  in 
Equation  3.4 


=  2A  +  2Bx  +  3C(x2  -  y2)  +  4D{x3  -  3xy2 ) 
6x 

=  A  -  2By  -  6Cx^  +  4D(y3  -  3x2y) 

By  Cauchy-Riemann  relationships 

<54  _  _64 

<5y  _  5x 

6 4  _  <54 

6x  ~~6y 

Equation  3.4  is  easily  evaluated. 


(44) 


3.3  Areas  of  Cells 


The  Jacobian  transformation  provides  a  convenient  method 
for  finding  the  areas  of  individual  cells  of  a  flow 
net.  Alternative  methods  are  based  on  co-ordinate 
geometry  via  subdivision  and  Gaussian  quadrature. 

Because  we  require  an  accurate  assessment  of  the  area 
to  calculate  the  quantities  appearing  in  Equations 
2.10  and  2.11  an  investigation  of  the  accuracy  of 
the  various  methods  for  finding  area  was  undertaken. 

The  flow  net  shown  in  Fig.  3.4  was  adopted  for  this 
exercise . 


Fig.  3.4  Orthogonal  net  for  Cell  Area 
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The  Jacobian  method  follows  from  the  expression  for  the 
the  area: 

A  =  |j  |*  Ax  *  Ay 

where  J  is  the  Jacobian  taking  on  the  lowest  values 
of  the  (x,y)  co-ordinates  of  the  given  cell.  The 
maximum  values  of  Ax  (and  Ay)  equals  unity  so 
convergence  can  be  tested  by  subdividing  the  cell  into 
smaller  areas,  say  Ax  =  0.1,  Ay  =  o.l,  and  summing  the 
subareas. 

Gaussian  quadrature  using  two  and  three  point  inter¬ 
polation  is  the  most  commonly  used  numerical  scheme. 

The  weighting  factors  are  given  in  Table  III-2  and  the 
area  is  obtained  from  _he  expression 

m  m 

A  =  £  I  w  *w.*  F (a . ,b . )  +  E  3.6 

j=l  i=l  1  3  13 

where  m  =  order  of  the  integration  rule 
w^  =  weights 

b^  =  abscissae  of  integration  points 
F  =  function  values 


E  =  error 

in  approximation 

TABLE  III-2 

Abscissas 

and  weights  for  Legendre- 

Quadrature 

( of  order  4 )  . 

Interval  |-1, 

1| 

m 

Abscissae 

Weight 

2 

+0.557350 

1 

3 

0 

8/9 

+0.7774597 

5/9 

4 

+0.339981 

0.652145 

+0.861136 

0.347855 
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Interval  )o,+lj* 


0.06943 

0.34785 

0. 33000 

0.65214 

0.69999 

0.65214 

0.93056 

0. 34785 

0.04691 

0.23692 

0.23076 

0.47862 

0. 5000 

0.56888 

0.76933 

0.47862 

0.95308 

0.23692 

For  the  curvilinear  square  shown  in  Fig.  3.4  with  origin 
at  x  =  7,  y  =  5  the  four  point  quadrature  over  the 
interval  |o,  +1 j  gives  an  area  of  442.60  units. 


The  results  for  cubic  interpolation  using  Equation  3.5 
are  as  follows 


Number  of  Subdivisions 
Nil 
2 

4 

5 

10 


Area  of  Cell 
388.00 
415.00 
428.75 
431.52 
437.08 


By  taking  the  Jacobian  determinant  at  the  point  x  =  7.5, 
y  =  5.5  (and  no  subdivisions)  the  value  obtained  for  the 
area  is  442.00  units. 


*  The  Legendre-Gauss  coefficients  for  the  interval  (0,+l| 
were  calculated  by  Al-Salihi  (1978). 
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Hence  it  appears  from  the  quantities  underlined  in  above 
that  the  Jacobian  determinant  evaluated  at  the  intersection 
of  the  medians  to  the  sides  gives  a  good  estimate  for  the 
asymptotic  value  of  the  area  projected  by  subdividing  the 
cell.  Where  tho  Jacobian  can  be  determined  analytically 
it  proves  more  expeditious  than  quadrature  formulae  which 
in  itself  is  a  motivation  for  trying  to  fit  polynomials 
to  flow  nets. 

3.4  The  Conjugate  Function 

The  foregoing  theory  pertains  to  conformal  mapping  from 
one  plane  to  another.  At  this  juncture  it  is  relevant 
to  compare  the  geometrical  exercise  with  its  physical 
counterparts,  namely,  the  derivation  of  velocity  potential 
and  stream  functions  in  fluid  mechanics,  or  the  plotting 
of  equipotential  and  flow  lines  in  electrostatics  and 
soil  mechanics.  In  these  problems  it  often  transpires 
that  one  of  the  functions  is  more  readily  plotted  than 
the  other  on  physical  grounds.  For  instance  as  an 
example  of  a  potential  function  let  the  distribution  of 
the  potential  be  denoted  by  <p  (which  is  analogous  to  the 
previous  interpretation  of  |  as  a  curve  in  the  plane  of 
Cartesian  space)  where  <f>  is  a  known  function  of  position 
in  the  solution  domain  in  the  x-y  plane.  In  this 
instance  <J>  represents  the  distribution  of  a  potential 
such  as  pore  pressure,  voltage  or  magnetic  flux.  Take 
for  example  the  real  part  of  the  harmonic  function 

<p  =  x2  +  4x  -  y2  +  2y 
where  <$>  satisfied  Laplace's  equation. 


A ii  +  ill 

<5x2  6y2 


0 
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The  conjugate  function  must  satisfy 


+  i!±  = 

6x2  6y2 


The  function  ip  is  deduced  from  the  expression 


In  this  example 

ry  rx 

^  =  2  \  (x+2) dy  +21  | (y-1) |  dx 

Jy  Jx„  y  yo 

J  o  o 

=  4y  -  2x  +  2xy  +  c 

where  c  =  2(x  -xy  -  2y  ) 

o  o  o  o 

Thusly  the  real  and  imaginary  parts  of  an  analytic 
function  are  evaluated.  However,  this  simple  example 
merely  serves  to  demonstrate  the  method  of  conjugate 
functions;  in  practical  examples  the  process  can  become 
very  difficult.  A  general  method  using  finite  difference 
approximations  on  a  curvilinear  orthogonal  grid  combined 
with  an  analytic  solution  has  been  reported  by  Centurioni 
and  Viviani  (1975) . 

The  discource  in  this  Chapter  leads  to  the  conclusion 
that  mapping  functions  can  be  deduced  by  inspection  of 
the  constraints  in  a  problem.  The  initial  trials 
can  be  adjusted  to  approximate  flow  paths  in  the 
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areas  of  fluid  and  soil  mechanics,  diffusion  process 
and  electrostatics.  Experience  of  generating 
orthogonal  nets  by  simply  varying  the  parameters  in  a 
cubic  polynomial  representation  suggests  that  a  useful 
technique  has  been  introduced.  By  using  the  computer 
generated  plots  to  overlay  a  given  potential  field  the 
functional  relationships  can  be  determined  piecewise,  or 
in  favourable  circumstances  the  mathematical  description 
of  the  entire  field. 

We  have  performed  these  exercises  on  a  variety  of  flow 
nets  with  the  result  that  it  has  always  been  found 
possible  to  locate  the  origin  if  the  net  had  been 
accurately  drawn  from  a  single  origin.  The  values 
listed  in  Table  II1-I  indicate  that  the  sets  are  distinct, 
leading  to  well  conditioned  simultaneous  equations.  The 
main  drawback  to  this  approach  is  the  stipulation  of  a 
unique  origin  of  the  net.  It  will  be  seen  in  the 
further  development  of  our  analysis  that  this  requirement 
can  be  relaxed  albeit  at  the  expense  of  computational 
effort.  A  further  drawback  is  due  to  the  fact  that  the 
conjugate  function  V  is  not  an  even  function,  as  may  be 
seen  by  inspection  of  Equation  3.1.  Therefore  it  results 
in  plots  with  skew  symmetry  which  limits  class  of  problems 
to  which  the  nets  may  be  applied.  The  only  symmetric  plot 
results  from  the  quadratic  terms  of  Equation  3.1. 

Negative  integers  employed  as  the  exponents  in  the  analytic 
function  lead  to  cumbersome  algebra  in  the  Jacobian  determinant, 
and  hence  taking  n  <  0  confers  no  practical  advantage. 
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CHAPTER  4 


Cubic  Interpolatory  Splines 


4.1  Introduction 


The  theory  of  splines  is  well  documented  so  only  a  brief 
review  iS  presented  herein  as  a  background  to  the  develop¬ 
ment  of  a  computer  program  for  parametric  cubic  spline 
interpolation  of  material  flow  paths. 

Splines  are  an  important  tool  in  modern  numerical  analysis 
for 

1)  Interpolation 

2)  Numerical  integration  and  differentiation 

3)  Numerical  solutions  of  ordinary  differential 
equations 

4)  Numerical  solutions  of  partial  differential 
equations. 

However  the  interpolatory  property  of  the  splines  is  only 
of  interest  here.  The  great  advantage  of  splines  in 
interpolation  is  that  they  do  not  have  the  oscillatory 
property  of  the  interpolating  polynomial.  The  most 
widely  used  splines  are  the  cubic  splines. 

4.2  Definition  of  a  Cubic  Spline 

A  cubic  spline  S(t)  with  modes  t. ,  t,  .  t 

i  ^  n 

(t^  <  t2  <  t^  .  <  tn)  is  a  function  which  in  the 

interval  t^  £  t  <  t^  +  1  (i  *  1*2,  ...  n-1)  reduced  to  a 
cubic  polynomial  in  t,  and  at  each  interior  node  t^ 

(i  =  2, 2,... n-1),  S (t) ,  S'(t)  and  S''(t)  are  continuous 
function  and  its  derivatives. 
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Such  a  function  may  be  represented  in  the  interval 
ti  <  t  <  ti  +  1  (i  =  1,2,... n-1)  in  the  form 

S±  =  F(t)  =  ait3  +  b±t2  +  cit  +  d± 

where,  associated  with  each  cubic  arc  S^f  there  are  four 
unknown  coefficients  a^,  b^,  and  d^.  it  follows  that 
there  are  4  (n-1)  conditions  to  be  fulfilled,  so  that  S(t) 
is  fully  defined  mathematically 

1)  S (t)  must  match  at  the  interval  nodes 
S(ti)  =  T(t±) 

i.e.  n  equations 

2)  S (t)  must  be  continuous  over  the  boundaries 
SUjJ  =  S(ti+)  i  =  2  (1)  n-1 

i.e.  n-2  equations 

3)  S'(t)  and  S''(t)  must  be  continuous  over  the 
boundaries 

S 1 (t±_)  =  S' (ti+) 

S'  '  (t±_)  =  S"  (ti+)  i  =  2  (1) n-1 

i.e.  2n-4  equations 

Hence  there  are  4n-6  equations  with  4n-4  unknowns.  To 
overcome  this  problem  the  natural  cubic  spline  can  be 
employed  which  has  the  extra  condition  for  the  second 
derivatives 

4)  S' ' (t^)  =  0 

s"(tn)  =  0 

This  corresponds  to  letting  the  physical  spline  project 
past  the  end  weights.  Usually  one  is  not  interested  in 
t  <  t^  or  t  >  t  but  if  they  arise  the  natural  spline  is 
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r 


r 


extended  by  straight  lines  agreeing  in  value  and  slope 

at  the  ends  t.  and  t  . 

i  n 

The  second  derivative  is  then  also  continuous  at  the  ends. 
Hence  there  are  two  more  equations  to  deduce  4n-4  equations 
to  solve  for  4n-4  unknowns.  Thus  the  mathematical  spline 
is  explicitely  defined. 

Condition  (1)  implies  that  the  curve  goes  through  the  knots 
or  nodes.  Conditions  (2)  and  (3)  imply  continuity  of 
alignment/  continuity  of  slope  and  continuity  of 
curvature  respectively. 

In  contrast  to  polynomial  interpolation  which  increases  the 
degree  to  interpolate  more  points,  here  the  degree  is  fixed 
and  one  uses  more  polynomials  instead  -  one  for  each 
interval.  For  each  interval  the  natural  cubic  spline  is 
unique  and  is  the  smoothest  interpolating  curve  through 
these  points. 


4.3  Derivation  of  the  Natural  Cubic  Spline 


For  notational  convenience  let 


=  fci+l  "  ti 


S1(t)  =  S[Z) ,  te  t±,ti+L 


=  S'*^) 

Because  S^(t)  is  a  cubic  polynomial  the  second  derivative 
S! '  (t)  is  a  linear  polynomial  and  can  be  expressed  in  the 


FIGURE  4.1 


M  M 

si,(t)  -  K7  ‘'irt-1'  +  irf  't-V  (4-l> 

i  =  1 (l)n-l 

To  obtain  expressions  for  S^(t)  equation  (4.1)  is 
integrated  twice  giving 

si(t>  -  SKj"  +  +  «  +  B 

where  A  and  B  are  constants 


To  determine  A  and  B  conditions  (1)  states  that 


Si(t1)  «=  F(ti) 


Si(titl> 


M 


Si(ti)  =  F(ti)  a  ^i+l'^i 


F(ti+i) 

i  ^  ,  Mi+l(ti“ti)* 


ltJAl-t,)3  + 


6h. 


+  Att  +  B 


Mihi2 

Si(ti)  =  F(ti)  =  ~-g-—  +  At±  +  B 


M  M 

si<W  -  P(ti+1>  *  irVm''  +  -arr'WV  ’  +  Rti+i+B 


154) 


or 


2 


By  linear  interpolation  (Fig.  4.2) 

M . h . 2  (t.,,-t)  Mi+lhi2 

At  +  B  «  (Fltj)  -  — ij-i -  +  <mi+1>  -  - ) 

(t-ti) 


This  leads  to  the  following  expression  for  S^t) 

M.  M. M.h .  ( t .+ . _t) 

sl(t)  -  +  +  -5— >— K7 - 


Mi+lhi  (t— t. ) 

+  <F«W  sf- 


(4.2) 


i  =  1,2,  ..  (n-1) 


The  condition  of  slope  continuity  leads  to  a  recursive 
equation  for  the  unknowns  i.e.  the  second  derivatives. 
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fl 


To  find  M,  .  the  condition  that  the  first  derivative  of 
( s ) 

S ^ ( t )  is  continuous  throughout  is  used 


i.e.  S  * i_1 (ti)  =  S*i(ti) 


Differentiating  equation  (4.2)  yields 
M, 

FF_ . 

i 


i  =  2 , 3 . .n 


F(ti+i,‘F(ti) 


s'i(tl  *  7K7<’Wt>2  +  +  - K. 


-  T<Mi+rMi> 


S' 


M,  .  M.  F (t.) — F (t._1) 

i-1  (ti)  =  2FiJ1’^ti"ti^2  +  2hi_1^ti"ti-l^  +  h^~[ 

-  — 

F<ti)-F,ti-i>  hi-i, 


M. 

_ i _ h2  + 

2hi-IH 


li-l 


-g — (Mi-Mi-l)  (4*3) 


M,  F ( t . , . ) -F ( t . )  h. 

s,i(ti}  =  ‘  Ih^(ti+i'ti) 2  +  - - g^itl'”^ 


M.  F(t.  ..)-F(tJ  h^ 

—Lh2  +  _  1+1 

2hiI  i  .»i 


i ) — F ( t . }  h. 

tt — -  '  -r(Mi+rMi) 


(4.4) 


Let  F(t±)  =  Fi 


Equating  (4.3)  to  (4.4)  gives 


M.  h 


F,-F4 


i“i-l  .  ‘i  xi-l  "i-1,..  ..  ^ 

+  “h^ - T-^i  Mi-1> 


Mih;  f.  ,-f.  hi 

-  - 


which  gives 
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<hi-l> 


M 


i-l  FT 


+  2M 


(hi+hi_i) 


+  M 


i+1 


F  -F 

r  i+l  i 


hi-l  1 


i  =  2,3  ..  (n-1) 


By  definition  of  the  natural  cubic  spline  and  Mn  are 
zero.  For  the  intervals  and  (tn_^,t  )  the  above 

equation  reduces  to 


2M. 


(h2+h1) 

“sr 


+ 


F  -F 
1 2  *1< 


M 


(hn-2j 

1  1 
!  +  2Mn-l| 

!hn-l'hn-2) 

1  6  1 

l  -  b  1 

[Vpn-1 

>h  .  < 

>  h  < 

h  ' 

>  h  *1 

(  n-1) 

1  ( 

{  n-1  , 

1  n-1  ( 

[  n-1 

F  -F 
n-1  rn-2< 


n-2 


Writing  the  equation  for  all  the  intervals  gives  the 
Equations  4.5  as  shown  on  the  next  page. 
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This  gives  a  system  of  tri-diagonal  linear  equations 

which  are  diagonally  dominent.  This  system  of  equations 

may  be  solved  by  using  the  Thomas  alogrithm  or  by  any 

method  for  solving  linear  equations.  However  one  can 

generate  a  recursive  algorithm  to  solve  for  the  M  .  Sampine 

s 

and  Allen  (1973) . 


Let  d . 
x 


6  (Fi+l~Fi  Fi'Fi-l) 
hi  hi  '  hi-l 


and  let  M.  .  =  p.M.  +  t.  i  =  2,3,...  n 

l-l  l  i  i 


(4.6) 


Whero  the  values  for  p^  and  must  be  found. 


=  O  therefore  one  may  assume  that  p^  and  =  0. 


Substituting  (4.6)  into  (4.5)  gives 


+  2[1+  ^7 lvMi+i  *  di 


or  {-E-ipi  *  2  1+-F^  >  Mi+Ml+1  -  di*  -KTTi 
1  (  1  )  i 


Therefore 


M.  = 

l 


-Mi+1 


h.  , 
d - 

h .  l 
l 


li-l.  .  J,,hi-l}  hi-l.  .  J,,hi-1 
“ pi  +  2  ~  1  mT- 


This  has  the  same  form  as  equation  (4.6)  hence 


pi+l 


i-l  .  .  i-l< 

— pi  +  21  +  — < 


(4.7) 
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and 


If  it  can  be  proved  that  no  denominator  vanishes  this 
gives  a  simple  recursion  for  computing  and 
(i  =  2,3,...,n).  Then  starting  with  Mn  =  0  equation  (4.6) 
may  be  used  to  calculate  in  the  order  (n-1,  n-2,  . . . ,2)  . 

hi-1  hi-l 

To  prove  — r-  --p .  +  2  (1h — r — )  ?  0  for  all  i 
ni  1  ni 

Clearly  I p^ I  <  1  If  I  pj  <  1  then 


2 

i 


>  2. 


Therefore  equation  (4.7)  holds  and  the  recursion  expression 
for  the  Mg  is  finite. 


4.4  Parametria  Cubic  Spline 

It  has  been  shown  that  for  a  variable  t  and  its 
corresponding  dependant  variable  F(t)  that  a  smooth  curve 
may  be  generated  to  interpolate  the  fixed  points. 

If  the  variable  t  is  assumed  to  be  the  variable  x  and  F(t) 
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assumed  to  be  y  (x  and  y  being  the  cartesian  plane  co¬ 
ordinates)  ,  then  for  any  x  value  using  Splines  S^(x)  may 
be  found.  But  this  is  of  no  direct  avail  in  representing 
a  multi-valued  function  of  x,  eg.  a  closed  curve,  or  a 
curve  that  doubles  back  on  itself  or  has  double  or  multiple 
points  or  if  it  is  an  open  curve  with  very  large  slope, 
dy/dx->-°°  at  some  point  in  its  range .  However  such 
functions  y  =  F(x)  can  be  represented  in  parametric  form. 


y  =  F  (t)  ;  x  =  G  (t) 


where  t  is  a  parameter  with  values  in  a  certain  interval; 
hence  the  two-dimensional  parametric  natural  cubic  spline 
is 


Fi(t) 
G±  (t) 


y  =  ait3+bit2+cit+di  (a) 

x  =  3+fit  3+git+hi  (b) 


(4.8) 


This  produces  two  splines  in  t  and  the  previous  theory  of 
splines  can  be  applied  to  each,  considering  nodes  in  the 
interval  of  t  which  are  strictly  increasing  in  value. 

Here  dy/dx  = 

Values  of  t  for  which  dx/dt  =  0  correspond  to  where  the 
slope  is  infinite.  In  this  type  of  cubic  spline  one  can 
find  in  general,  three  separate  values  of  t  which  give 
the  same  value  of  x  and  of  course  in  general  three 
different  y  values  allowing  it  to  represent  curves  which 
form  loops  or  even  double  loops.  It  has  been  specified 
that  t  must  be  strictly  increasing  between  (x^y^  and 
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xn»yn) »  that  is,  on  the  curve  in  the  x,  y  plane  if  t  takes 
the  value  at  the  point  (x^,y^)  and  taken  the  value  t 
at  the  point  (x,y)  and  takes  the  value  T^+^  at  the  point 
^xi+l'^i+l^  t^ien  Ti  <  Ti+i*  This  means  that  in  the  (t,x) 
plane  and  (t,y)  plane  x  and  y  are  single  valued  functions 
of  t.  If  this  parameter  t  is  considered  as  the  distance 
along  the  curve  then  the  resulting  cubics  in  t,  for  x  and 
y  give  the  curve  in  the  (x,y)  plane  of  that  formed  by  the 
spline.  If  t  is  the  distance  from  the  point  (x^,y^)  to 
the  point  (x,y)  on  the  spline  along  the  spline  then  the 
spline  may  be  represented  by  equations  (4.8)  (a)  and  (b) . 
The  values  T^  at  the  point  (x^y^  are  not  known  and  must 
be  estimated. 

4.4.1  Estimation  of  Curve  longth 

As  the  parameter  t  is  considered  the  distance  along  the 
curve,  then  the  curve  length  must  be  estimated  before  the 
curve  is  found.  The  procedure  used  is  to  write  a  cubic 
polynomial  of  the  form 

t  =  ax3  +  bx2  +  cx  +  d  (4.9) 

for  each  interval  |x^,  |  and  evaluate  the  integral 

L  =  Jx *i+1  {1  +  (dy/dx) 2 }^dx 

The  integration  is  formed  numerically  using  Gaussian 
quadrature  with  four  Gauss  points. 

The  curve  must  be  well  behaved  that  is  9^  or  9^+^  must  not 
equal  90°  and  10^  -  <  180°;  if  these  two  conditions 


(62) 


are  not  satisfied  the  line  must  be  rotated  towards  the 
horizontal  before  integration. 

The  resulting  lengths  for  each  interval  are  added  to  give 
the  length  of  the  whole  curve.  Obviously  as  this  is  an 
estimated  length,  the  error  accumulates  as  the  curve  length 
increases.  When  the  spline  is  found  using  these  estimated 
lengths,  the  lengths  are  re-calculated  and  by  an  iterative 
process  an  improved  curve  length  may  be  found  hence  giving 
a  better  fit  to  the  locus  of  the  knots. 

4.4.2  Estimation  of  Slope  at  node  points 

An  alogrithm  to  provide  values  for  the  slope  at  the  point 
can  be  developed  by  finding  the  equation  of  a  curve 
through  P^^  and  several  adjacent  points  and  by  taking  the 
derivative  of  the  curve  at  P^,  McConologue  (1970) .  The 
slope  at  each  node  is  needed  to  determine  the  constants 
a,  b,  c  and  d  in  equation  (4.9). 

To  find  these  four  unknowns,  four  equations  are  needed 
within  the  interval  (x^  xi+1)  . 

These  are 

^i  =  axi  +  kxi  +  cxi  +  ^  (4.10) 

yi+l  =  axi+l  +  bxi+l  +  CXi+l  +  d 
tan6^  =  3ax*  +  2bx^  +  c 

tan©i+1  =  3ax£+1  +  2bxi+1  +  c 
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The  function  to  be  differentiated  is  in  parametric  form 
to  be  compatible  with  the  parametric  spline.  The  method 
used  is  to  pass  a  parabola  through  three  successive  points 
Pi-lf  p^,  Pi+^»  the  x  and  y  co-ordinates  being  given 
independently  in  terms  of  a  parameter  u.  Such  a  parabola 
is  not  unique,  but  values  for  the  slope  acceptable  over  a 
wide  range  of  applications  is  obtained  by  writing  x  and  y 
independently  as  Lagrangian  interpolating  polynomials  in 
u,  where  x  and  y  have  the  values  and  yA_^  at 


u  =  Di-1'  xi'*i 

at 

u  =  0  and  xi+1,  yi+1  at  u  = 

where 

Di 

= 

{(Ax.^2  +  (Ayi)2}!|s 

(4.11) 

Ax^ 

= 

xi+l  ~  xi 

Ayi 

= 

*i+l  ‘  yi 

Differentiating  these  equations  and  simplifying  gives  the 
following  equations  for  the  slope  at  any  point  pr: 


where 


Cos6r 

=  Cr/Nr 

Sin6r 

■  Sr/Nr 

(4.12) 

Nr  " 
slope  = 


Axi-lar  +  Axier 


Ayi-lar  +  Ayi6r 

c2  +  S2  ^ 
r  r 


(4.13) 
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On  the  assumption  that  the  positive  direction  along  the 
curve  is  from  to  Pi+1* 


and 


a 


r 

r 


r 

r 


Di^2Di-l  +  Di )>  D2i  and  “°2i  for 
i-1/  i,  and  i+1  respectively 


i-1 ' 


i-1*  Di-l^Di-l 


+  2D.  ) 


for 


i-1,  i  and  i+1  respectively. 


The  form  r  =  i  is  normally  used  where  possible,  those  for 
r  =  i-1  and  r  =  i+1  being  used  for  the  beginning  and  end 
of  an  open  curve  or  at  points  on  an  open  or  closed  curve 
where  there  are  cusps  and  other  discontinuities  in  the 
derivative.  However  the  method  used  here  in  obtaining 
slopes  at  node  1  and  node  n  is  to  generate  nodes  0  and 
n+1  by  taking  xq  as  x^  -  (x^x^)  and  xr+1  as  xn  +  (xn“xn_^^ 
and  by  spline  interpolation  finding  the  corresponding  y 
values  i.e.  yQ  and  yn+1-  Hence  nodes  1  and  n  are 
internal  nodes  thus  the  form  r  =  i  is  used. 


4.4.3  Interpolating  Intermediate  Points 

At  the  nodes  the  value  of  the  parameter  t  coincides  with 
the  distance  from  an  arbitrary  point  (x^,y^) .  It  is  not 
true  that  an  intermediate  value  of  the  parameter  t 
represents  a  point  that  has  the  same  value  for  the 
distance  to  the  point,  for  if  distance  {(x^,y^),  (x,y) } 
equals  the  distance  along  spline  from  (x^,y^)  to  (x,y) 
and  if  distance  {(x^,y1),  (x^jy^)  }  =  s.^  and  is 

a  node  point  with  the  parameter  value  t  =  t^ 


JtMa  — 1 


(6b) 


V 


then  t^  = 

But  if  distance  {(x^y^).  (x,y) }  is  s  and  (x,y)  is  some 
point  on  the  spline  with  parameter  value  t  *  T  j*  t^  then 
T  =  s  is  not  necessarily  true.  However  there  is 
necessarily  continuity  of  the  parameter  t  and  distance  s 
and  the  parameter  values  are  coincident  with  the  distance 
values  at  the  nodes.  From  this  it  is  known  between  which 
two  nodes  a  given  value  of  the  distance  s  lies,  the  bounds 
between  which  the  associated  parameter  value  T  will  lie 
can  be  found. 


If  si  <  s  <  si+l  then  t  <  T  c  t1+1 

In  terms  of  the  parameter  t  in  equations  (4.8)  (a)  and 
(b)  the  distance  along  the  spline  is  given  by 

s  =  t±  +  (F‘(t)}2  +  (G[(t)2}*dt  (4.14) 

The  numerical  evaluation  of  this  expression  for  the 
unknown  T  yields  the  associated  parameter  of  the  distance 
s,  and  hence  x  and  y  coordinates  for  any  point  on  the 
curve . 

The  procedure  is  illustrated  schematically  in  the  mapping 
shown  in  Fig.  4.3. 
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FIGURE  4.3:  Parametric  Mapping  of  plane  curve 


4.5  Data  Input 

A  computer  program  based  on  foregoing  theory  was  written 
in  Fortran  IV  for  interpolating  flow  nets  by  parametric 
cubic  splines  (O'Laoide,  1979).  The  computer  program 
is  structured  to  compute  x  and  y  co-ordinates  for  given 
lengths  along  a  curve.  It  also  computes  x  co-ordinates 
for  given  y  co-ordinates  and  calculates  the  slope  of 
the  curve  at  a  given  x  co-ordinate. 

The  number  of  nodes  with  their  co-ordinates  must  be  read 
inland  when  the  number  of  nodes  is  given  as  zero  the 
program  is  terminated. 

The  program  has  six  input  types  and  are  as  follows: 

Two  pointers,  denoted  by  L  and  M,  enable  selection  of  a 
set  of  options  as  follows: 
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L=1 


a  length  is  read  in  and  the  corresponding  x  and 
y  co-ordinates  are  calculated 


L=2  a  length  of  increment  is  read  in  with  an  upper  and 
lower  bound.  The  x  and  y  co-ordinates  for  each 
incremental  length  are  calculated 

L=3  an  x  co-ordinate  is  read  in  and  its  corresponding 
y  co-ordinate  is  calculated 

L=4  an  x  co-ordinate  increment  is  read  in  with  an  upper 
and  lower  bound.  The  y  co-ordinates  for  each 
incremental  x  co-ordinate  are  calculated 

L=5  an  x  co-ordinate  is  read  in  and  the  slope  of  the 
curve  at  that  co-ordinate  is  calculated 

L=6  an  x  co-ordinate  increment  is  read  in  with  an  upper 
and  lower  bound,  and  the  corresponding  slopes 
within  the  range  are  calculated. 

When  M  =  0,  for  L  =  1,2,3  and  4  the  slopes  corresponding 
to  the  x  co-ordinates  will  also  be  calculated.  When  M  =  2 
a  cumulative  chord  length  is  read  in  and  the  corresponding 
x  and  y  coordinates  are  calculated. 

4.6  Applications  of  Cubic  Spline  in  Harmonic  Analysis 
Interpolation  of  Contour  Intervals 

Normally  a  spline  is  viewed  as  an  ordinary  function  which 
gives  the  value  of  the  dependent  variable  for  a  pre¬ 
scribed  value  of  the  independent  variable  (in  a  closed 
interval).  In  what  follows  it  is  demonstrated  that  a 
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reverse  process  can  be  implemented  whereby  the  dependent 
variable  can  take  on  prescribed  values  and  the  original 
independent  variable  can  be  treated  as  a  dependent 
variable.  This  transition  provides  a  means  of  plotting 
contours  for  monotonic  functions  and  in  particular 
solutions  of  Lapace's  equation  with  coordinates  (x,y, 

0)  say,  where  0  is  the  potential. 

To  map  the  contours  we  require  ordinates  at  stations 

along  a  grid  line  where  the  ordinates  differ  by  a  preset 

value  i.e.  the  contour  interval.  A  cross-section  on 

any  grid  line  will  normally  have  contour  ordinates  at 

irregularly  spaced  distances  from  the  origin  of  the 

cross  section.  Each  contour  ordinate  t  is  located 

c 

at  a  root  of  the  equation 

3.  U)  -  tc  =  0  4.15 

wnere  is  the  spline  function  representing  f(t) 

through  the  knots  in  the  interval  jx^,  x.^  +  l|  or  in 
interval  | y ^ ,  yi  +  lj. 

The  roots  of  the  Equation  4.15  can  be  determined  with 
any  degree  of  accuracy  by  repeated  application  of  the 
Newton-Raphson  formula,  provided  that  there  is  no  local 
extremum  of  the  function  in  the  neighbourhood  of  a  root. 
This  condition  is  satisfied  by  all  harmonic  functions 
if  sources  and  sinks  are  isolated  in  the  physical 
problem. 

Reverting  to  the  notation  of  the  parametric  cubic  spline 
Equation  4.2  and  4.15  are  written  in  the  form: 
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fci+l  -  fci 


4.16 


3y  applying  the  Newton-Raphson  formula 


S.  It) 

t  =  t  -  _ - _ 2 _ 

n+1  n  S[(t  )  4.17 

n  =  iteration  number 

the  first  approximation  to  a  root  within  the  interval 
| t ^ ,  ti  +  jJ  is  conveniently  derived  by  letting  t^  =  t^, 
or  alternatively  t  =  t  . ..  The  roots  tor  a  given 
value  tc  will  be  available  on  lines  parallel  to  both 
x  and  y  coordinate  axis  by  simple  looping.  However  it 
is  necessary  to  sort  the  array  of  roots  to  ensure  that 
those  roots  in  close  proximity  are  grouped  together  for 
machine  plotting  of  the  contour  ordinate  as  a  sequential 
curve . 


Transformation  of  Boundary  Conditions 

To  fully  avail  of  the  properties  of  harmonic  functions 
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an  automatic  means  of  converting  boundary  conditions  is 
required.  The  Dirichlet  problem  is  one  or  finding  a 
function  u  to  satisfy 

y2u  =  0  (u)  =  f  4.18 

where  f  is  a  function  of  position  on  the  boundary  s  and 
V 2  is  the  Laplacian  operator. 

The  Neumann  problem  is  that  of  finding  a  function  u  to 
satisfy 


where  g  is  a  given  function  of  position  on  the  boundary 

s  and  is  the  normal  derivative  outward  from  the 
Sn 

region  considered.  By  using  the  Cauchy-Riemann 
equations  either  the  Dirichlet  or  the  Neumann  boundary 
value  problems  may  be  transformed  into  the  other. 

On  the  boundary  s  the  Cauchy-Riemann  equations  imply 

6u  _  _  _6v 

6  s  6n 

where  yy  is  the  tangential  derivative  of  the  function  u 
in  the  positive  sense  (region  of  analyticity  on  the  left 
i.e.  counterclockwise  circulation)  along  the  boundary, 
and  yy  is  the  normal  derivative  of  the  conjugate 
function  outward  from  the  region.  From  analyticity 
requirements  the  integral  of  yy  around  a  complete 
boundary  must  be  zero  (i.e.  property  of  perfect 
differential;  so  the  Neumann  problem  (yy)  =  g  has  no 
solution  unless  g  satisfies 
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o 


4.20 


£ s9  ds  = 


The  cubic  spline  provides  a  means  of  checking  boundary 
conditions  via  the  first  derivatives  at  tne  knots. 


4.7  Hilt 'RES ENT  A T  ION  OF  SPLINE  FUNCTIONS  BY  TRUNCATED  POWER  SERIES 


2 

The  cubic  spline  enforces  C  -  continuity  at  the  knots  and 
hence  gives  an  adequate  representation  of  the  mean  stress 
contours  of  the  harmonic  problem.  On  the  other  hand  the 

4 

solution  of  the  biharmonic  equation  demands  C  -  continuity, 
therefore  for  completeness  reference  is  made  to  the  truncated 
power  series  representation  known  as  the  Schoenberg-Whitney 
formula  as  follows: 


S (x )  =  P (x)  +  J  C  (x  -  x.)™ 
J=1  : 

where  P  e  P 

and  (x  -  x.)m  =  0  If  x  $  xT 

3  +  J 

=  x  -  x.  if  x  >  x. 

3  3 


4.21 


P  denotes  the  class  of  polynomials  of  degree  m  or  less,  and 

c - c  are  the  unknown  constants. 

1  n 

In  beharmonic  analysis  the  appropriate  spline  is: 

2345 

S (x )  =  a^  +  a. x  +  a0x  +  a^x  +  a.x  +  a^x  + 
o  1  2  3  4  5 

n  c; 

J  C  (x  -  x  ) 5  4.22 

J  =  1  J  J  + 

where  the  interval  is  divided  into  (n+1)  equal  length 

segments.  There  are  (n  +  6)  unknown  constants  in  the 

quintic  spline  which  are  evaluated  by  collocation  at  the 

n  knots  (x. ,  x_ - x  )  and  at  the  end  points  x  and 

12  n  o 

x^  +  ^  by  using  two  boundary  conditions  at  each  end  of  the 
interval.  The  Schoenberg-Whitney  formula  produces  a  cubic 
spline  if  m  =  3  and  the  set  of  polynomials  P  is  equal  to 
or  less  than  P  .  The  resulting  representation  is  an 
aggregate  of  B-splines. 
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CHAPTER  5 


Harmonic  Analysis  of  Stress  Fields 


5.1  Finite  Difference  Approximations 

As  mentioned  earlier  the  solution  of  Laplace's  Equation 
gives  the  distribution  of  mean  stress  within  a  two 
dimensional  domain  provided  that  precise  boundary 
conditions  can  be  enforced.  In  this  Chapter  the 
question  of  extending  harmonic  analysis  to  axisymmetric 
problems  is  investigated.  The  objective  is  to 
prescribe  sufficient  ooundary  conditions  to  enable  a 
numerical  solution.  The  choice  of  finite  difference, 
finite  element  or  boundary  element  techniques  is  one 
of  personal  preference;  the  three  methods  are  equally 
applicable  in  the  majority  of  cases  but  the  FEM  has 
some  special  advantages  for  non-homogeneous  problems, 
Fenner  (1975)  Frind  (1977)  Huebner  (1975)  and  several 
other  authors.  In  stress  analysis  these  methods  yield 
data  for  plotting  rhe  contours  of  mean  stress  if 
Dir ichlet-type  conditions  are  imposed  on  the  boundaries 
(Neumann-type  conditions  may  be  more  easily  specified 
in  other  areas  suen  as  fluid  mechanics) .  Where  the 
analysis  is  confined  to  a  subdomain  of  a  continuum  as 
in  the  axisymmetric  flow  problem  the  specification  of 
boundary  conditions  is  identical  tor  any  one  of  tne 
numerical  methods  (finite  differences,  fininte  elements 
or  boundary  elements) .  The  values  of  the  mean  normal 
stress  must  be  prescribed  at  the  external  nodes  of  the 
discretised  domain.  The  finite  difference  method  is 
adopted  herein  because  it  simplifies  the  numerical 
calculations  if  the  domain  is  assumed  to  consist  of  a 
homogeneous  and  isotropic  half  space  (where  it  is 
possible  to  specify  a  constant  grid  interval  on  mutually 
perpendicular  coordinate  axis) . 
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The  governing  dif terential  equations  for  axially 
symmetric  stress  distribution  take  the  form,  TimosnenKO 
and  Goodier  (1951); 
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where  0  denotes  the  mean  normal  stress  multiplied  by 
a  factor  of  three,  and  the  remaining  terms  are  in  the 
usual  notation  for  cylindrical  coordinates.  By 
adding  (a)  (b)  and  (c)  of  above  set  it  can  be 

shown  that  the  governing  equations  reduce  to  the 
following  expression: 


620  1  60 

or’  *  61’ 


5.2 


where  0  can  be  interpreted  as  a  stress  potential  in 
the  absence  of  body  forces. 


5.2  Boundary  Conditions  in  Terms  of  Mean  Stresses 

Dirichiet  Boundary  Value  Problem  The  stress 
distribution  in  the  elasto-isotropic  medium  for  a 
uniformly  loaded  circular  area  can  be  calculated  by 
means  of  the  Boussinesq  solution  for  a  concentrated  load, 
the  principle  of  superposition  and  Maxwell's  reciprocal 
theorem.  The  stresses  in  the  vertical  and  radial 
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f 


directions  are  given  by  the  expressions;  Tiirosheuko  and 
Goodier  (1951) : 
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where  P  is  a  concentrated  vertical  load  and  r,z,9  are 
the  ordinates  of  cylindrical  coordinate  system. 


The  stress  components  can  be  found  directly  but  the 
calculations  are  cumbersome,  Harr  (196b)  JumiKis  Il9b9) 
Ahlvin  and  Ulery  (1962) .  The  direct  methods  generally 
reduce  to  the  evaluation  of  elliptic  intergals  which 
are  not  easily  handled  in  a  computer  program.  In  this 
section  an  approximate  method  based  on  superposition 
is  presented. 


The  uniform  loading,  or  other  continuous  distribution, 
is  replaced  by  equivalent  line  loads  on  annular  rings. 
The  method  involves  summing  the  stresses  generated  at 
the  boundary  nodes  by  point  loads  on  differential 
arc  lengths. 
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Let  P  denote  tne  line  load  per  unit  length  of  annular 
ring  as  shown  in  Fig.  5.1  and  let  W  denote  the  total 
load  on  the  annular  ring. 


Then  P  =  W/2Tir  where  r  is  the  mean  radius  of  ring  and 

3.  cl 

the  distance  to  the  boundary  of  finite  difference  mesh  is 
given  by 

0  2  ** 
d  =  (r2+h  -2rh  Cos6) 


5.1  Normal  Load  over  Circular  Area 

The  stresses  on  the  edges  of  any  diametral  plane  passing 
througn  the  axis  of  symmetry  are  given  by  the  expressions 
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Because  of  symmetry  the  superposition  of  stresses  is 
performed  in  the  range  0  <_  9  <_  it  with  small  increments 
of  0  corresponding  to  unit  length  of  arc  for  each 
annular  ring.  Nodes  on  the  axis  of  symmetry  are 
treated  separately  because  the  stresses  can  be  related 
to  total  load  W  on  the  rings.  The  stress  components 
on  the  axis  of  symmetry  for  uniform  pressure  on  circular 
contact  area  are  given  by  the  expressions 
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where  qQ  denotes  tne  uniform  loading  distribution  over 
the  area  of  the  circle  of  radius  a. 


5.3  Boundary  Tractions 

Discontinuities  i:i  surface  tractions  such  as  occur 
at  the  extremities  of  a  uniformly  distributed  load  are 
detrimental  in  a  finite  difference  analysis.  Smoothing 
by  equivalent  continuous  functions  of  the  load 
aistrioution  seems  imperative.  A  uniform  load 
distribution  over  part  of  tne  boundary  is  reasonably 
approximated  by  the  expression 

q  =  qQ  Sech2  £ 

0  <_  p  £  h 

oecause  the  definite  integral  of  the  function  converges 
within  a  short  distance  of  the  edge  of  the  uniformly 
loaded  area ^as  shown  in  the  following  taDulation).i . e . 
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TABLE  V-l:  Approximating  Function  for  Symmetrical 
Square  Wave 


p 

a 

Sech  (— ) 

cl 

Sech2  (— ; 

cl 

tanh  i-£) 

cl 

0.0 

1 . 0000 

1 . 0000 

0 . 0000 

0.  2 

0. 9803 

0.  9610 

0.1974 

0.  4 

0. 9250 

0. 8556 

u. 3799 

0.  6 

0.8435 

0.7115 

0.5370 

0.8 

0.7477 

0.5591 

0.6640 

1.0 

0.6480 

u.  4199 

0.7616 

1.2 

0. 5523 

0. 3050 

0.8336 

1.4 

0.4649 

0. 2161 

0.8853 

1.6 

0. 3880 

u. 1505 

0.9217 

1.8 

0. 3218 

0. 1036 

0.9468 

2.0 

0.2658 

0.0706 

0. 9640 

2.4 

0. 1800 

0.0324 

0.9837 

3.0 

0.0993 

0.0093 

0.9950 

Accelerated  convergence  to  a  uniform  load  distribution 
can  be  achieved  by  assuming  weighting  factors  in  the 
range  1  to  2;  one  such  factor  leads  to  the  expression; 

q  =  q  Sech2  /2  —  5.6 

n  ^o  a 

wiiere  qQ  is  tne  amplitude  of  the  uniformly  distributed 
load  and  p  denotes  distance  from  axis  of  symmetry.  The 
potential  function  on  free  boundary  is  8  =  -q. 

The  Dirichlet  conditions  of  the  original  problem  are  thusly 
replaced  by  equivalent  distributions  of  mean  normal  stresses 
on  the  boundaries  with  nodal  values  taking  on  the  ordinate 
values  of  this  distribution  via  a  superposition  technique. 


5.4  Finite  Difference  Operator 

Let  the  differential  Equation  5.2 
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be  represented  by  difference  equations 
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By  differentiating  the  numerator  and  denominator  (as 
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where  the  nodes  and  intervals  are  shown  in  Fia.  5.2 
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Fig.  5.2  Finite  Difference;  Mesh 


Neumann  Boundary  Value  Problem  The  Equation  5.2 

for  the  potential  function  in  terms  of  mean  stress  is 
anaiagous  to  the  equation  of  velocity  potential  in 
fluid  mechanics.  The  starting  point  for  its 
derivation  is  the  Continuity  Equation  of  an  in¬ 
compressible  fluid. 
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v/here  V  denotes  velocity  ar.d  $  velocity  potential. 


By  specifying  boundary  conditions  based  on  the  first 
derivative  of  the  solution  for  the  stream  lines  is 


obtained.  The  Neumann  problem  is  posed  by  taking  as  the 
Continuity  Equation  the  expression 
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where  the  boundary  conditions  are  those  discussed  in  Chapter  4. 
The  solution  of  the  Neumann  problem  can  be  extracted  to  within 
an  arbitrary  constant  throughout  the  domain.  There  are 
sufficient  data  for  contouring  lines  of  constant  i|>,  and  in  fact 
the  same  soLution  routine  can  be  used  for  finding  $  and  <|> 
simply  by  inserting  a  linking  subroutine  to  generate  the  tp 
values  at  the  nodes  of  the  boundary.  An  interesting  variant 
of  this  approach  has  been  reported  for  the  problem  of  seepage 
through  an  earth  dam,  Shaug  and  Bruch  (1976).  However,  for  the 
stress  distribution  problem  the  Neumann  conditions  pertain  to 
strain  energy  density  which  are  difficult  to  prescribe.  Hence 
an  alternative  method  is  required  for  finding  orthogonal 
trajectories  to  the  mean  stress  contours. 


The  problem  of  orthogonal  slip-line  fields  is  extensively 
covered  in  applications  of  plasticity  theory,  hence  as  a 
preliminary  to  the  present  author's  approach  to  potential  field 
problems  a  brief  outline  of  mechods  employed  in  plastic  analysis 
is  included  in  the  next;  section. 
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5.5  Slip-Line  Fields  in  Plasticity  Theory 

Method  of  Characteristics 

The  theory  of  the  ideal  plastic  solid  yields  differential 
equations  for  solution  of  the  geometrical  constraints  on  the 
slip  lines.  For  the  case  of  plane  strain  the  governing 
differential  equation  is  derived  on  the  following  assumptions : - 

1.  the  shear  stress  has  a  yield  value  k  which  is  induced 
at  every  point  on  a  slip  line, 

2.  the  equilibrium  equations  are  identical  to  those  of 
elasticity  theory, 

3.  the  direct  stress  normal  to  the  slip  line  is  the  mean 
stress  devoted  by 

4.  the  geometrical  constraints  on  the  slip  lines  can  be 
referred  to  Cartesian  coordinates  with  the  angle  subtended 
by  normal  to  the  plane  of  maximum  shear  and  the  x-axis  the 
identifying  parameter  (denoted  by  a)  as  depicted  in  Fig.  5.3. 

Consider  a  homogeneous  weightless  material  in  a  state  of 
limiting  equilibrium.  The  equilibrium  equations  are  given  by 
the  expressions,  Unksov  (1961):- 


6x  6y 


0 


0 


5.10 


and  the  stress  state  in  terms  of  o  and  k  is  as  follows 


a  +  k  Sin  2a 
a  -  k  Sin  2a 
-  k  Cos  2a 


5.11 
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By  substituting  the  expressions  5.11  into  5.10  the 
equilibrium  equations  in  terms  of  the  mean  stress  become: 


P  +  2k  (Cos  2a  +  Sin  2a  4s-)  =0 

Sx  6x  6y 

p-  -  2k  (Cos  2a  p  -  Sin  2a  =0  5.12 

6y  <5y  <5x 


On  eliminating  a,  by  differentiating  the  first  of  5.11  with 
respect  to  x  and  the  second  to  y  the  single  differential 
equation  is  obtained  by  subtraction  viz. 


-  +  2  Cot 


&  a  6  2  a 
6x6y  6y2 


6  a 
6y 


+ 


2  Cot 


(— )  2 


-  ( 


—  )2! 
6x 


=  0 


5.13 


By  the  theory  of  differential  equation  the  so-called 
characteristic  equation  for  5.12  is  given  by 

-  dy^  -  2  Cot  2a  dx  dy  +  dx2  =0  5.14 


The  roots  of  equation  5.13  are  both  real  i.e.  the  equation 
5.12  is  hyperbolic,  and  lead  to  the  solution  for  the 
characteristic  directions  as  follows: 


P 1 ) 

'dx' 


2 

Cot  2a  +  (Cot  2a  + 

Cos  2a  1 

Sin  2a  Sin  2a 


1) 


tana 


Cot  2a  -  (Cot  2a  +  1) 


h  _ 


Cota 
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It  follows  that  the  tangents  to  the  slip  lines  at  any  point 
form  angles  with  the  x-axis  which  are  given  by  the  equalitie 
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tan  a 
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Cot  a 


*£  » 
dx 


In  other  words,  the  characteristics  coincide  with  the  slip 
lines  and  form  an  orthogonal  system;  the  characteristics  cut 
a  free  surface  at  a  constant  angle  of  45°.  Thus  the 
geometrical  constraint  is  that  the  slip  lines  should  form 
an  orthogonal  mesh  within  the  solution  domain. 

Because  the  foregoing  equation  is  of  the  hyperbolic  type 
it  is  theoretically  feasible  to  propagate  the  stress  field 
solution  by  integration  along  the  characteristics,  Harr 
(1966)  pp  24  to  253.  However,  the  procedure  is  difficult 
to  implement  because  it  requires  much  physical  interpreta¬ 
tion  as  demonstrated  by  Sokolovsky  and  other  researchers. 

On  the  other  hand,  purely  geometrical  techniques,  such  as 
Prager's  graphical  method,  is  favoured  for  intricate 
problems. 


Matrix  Technique  for  Constructing  Slip  Lines 

A  matrix  analysis  for  constructing  slip-line  solutions  has 
been  reported  by  Dewhurst  and  Collins  (1973) .  The  procedure 
is  based  on  a  power  series  representation  of  the  solution  to 
the  governing  equations  first  used  by  Ewing  (1967) .  By 
starting  with  base  slip-lines  with  radii  of  curvature  Ro(a) 
and  Sq(6)  the  initial  curvatures  are  expanded  as  power 
series  in  the  angular  coordinate. 


V°) 


00  n  «  0n 

=  J'a  &n  ^  ,So(e)  =  ^  bn  nT 

n=o  n=o 
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By  introducing  Mikhlin  co-ordinates  (x,y)  it  can  be  shown 
that  the  relationship  between  the  coefficients  in  power 
series  expansions  of  these  co-ordinates  and  of  the  radii 


fcn+l  "  fcn-l 


where  the  terms  are  defined  in  the  original  paper.  The 
subsequent  matrix  formulation  leads  to  a  general  method 
for  'marching  out'  the  solution  of  plane  strain  rigid 
perfect-plasticity  problems.* 

Of  course  the  finite  element  method  has  made  an  impact  in 
plasticity  problems.  Zienkiewicz  and  Godbole  (1975). 

These  authors  compare  various  methods  for  treating 
viscous  incompressible  flow  with  special  reference  to 
Non-Newtonian  (Plastic)  fluids. 

Essentially  all  the  methods  aspire  to  constructing  the  a-8 
net  which  gives  a  graphical  picture  of  the  deformation  mode. 


* 


Program  listing  in  FORTRAN  is  appended  to  the  paper, 
Dewhurst  and  Collins  (1973) . 


5.6  Orthogonal  Trajectories  to  Countoure  of  Mean  Streeaee 
A  New  Numerical  Method. 

In  the  course  of  this  study  it  became  evident  that  a 
technique  more  versatile  than  that  presented  in  Chapter  3 
is  needed.  In  particular  if  the  special  attributes  of  para¬ 
metric  splines  are  to  be  exploited  for  plotting  trajectories 
of  mean  and  shear  stresses  a  solution  format  analogous  to 
the  matrix  method  for  plasticity  is  a  desirable  goal.  Apart 
from  applications  of  complex  variable  theory  and  trial  and 
error  sketching  of  flow  nets  no  techniques  emerged  from  an 
extensive  literature  search.  The  present  author  expended 
considerable  effort  to  simply  find  a  framework  for  the 
problem,  because  there  are  a  number  of  possible  formulations 
for  instance,  such  as  differential  forms,  optimization 
theory  or  indeed  eigenvalue  solutions.  It  transpires  that 
optimization  theory  affords  the  most  tangible  approach  and 
is  the  basis  of  the  method  proposed  herein. 

Inspection  of  the  plots  of  orthogonal  nets  produced  by  the 
harmonic  analysis  of  Chapter  3  indicated  that  each  set  of 
trajectories  of  the  a- 6  net  are  the  envelopes  of  a  family 
of  circles  as  illustrated  in  Fig.  5.4.  By  assuming  that 
points  of  tangency  of  certain  members  of  the  infinite  set 
of  circles  bounded  by  the  mean  stress  contours  lie  on  the 
required  trajectory,  the  optimization  problem  for  finding 
the  locus  of  the  centres  reduces  to  the  statement: 

'Given  two  non  intersecting  curves  in  a  plane  find  the 
family  of  circles  that  meets  the  constraint  of  tangency 
to  the  given  curves  and  in  addition  ensures  that  no  one 
member  of  the  family  intersects  any  other  member. ' 
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In  other  words  the  optimization  problem  in  algebraic  terms 
is  as  follows: 


Minimize 

C 

=  (xx  -  hr  + 

Subject  to 

the 

equalities 

fl(xl) 

= 

Pi  (x^ ’^2} 

f 2 (x^) 

= 

P2 ^X1 ' x2 ^ 

f3<Xl> 

S 

P3  (x1,x2) 

Fig. 

where  the 

solution  domain  is 

Xt 


Posed  in  this  fashion  the  optimization  problem  is  nonlinear 
and  moreover  the  constraints  are  not  known  a  priori;  in 
particular  the  co-ordinates  of  the  constraint  represented 
by  P_(x.,x2)  depends  on  a  previous  solution  step.  One  way 
of  making  the  problem  tractable  is  to  consider  one  circle 
at  a  time  and  to  use  a  local  coordinate  system  suitably 
located  on  one  of  the  fixed  constraints.  An  algorithm  to 
establish  the  centre  and  minimum  radius  which  satisfies  the 
constraints  is  the  initial  requirement.  When  the  family  of 
circles  are  defined  by  the  co-ordinates  of  the  points  of 
tangency  the  required  trajectory  passes  through  these  points. 
In  a  numerical  solution  it  will  not  be  necessary  to  graph 
the  circles.  This  so-called  a- 3  net  is  orthogonal  but 
certainly  is  not  identical  to  the  a- 3  net  of  plastic 
analysis. 

Consider  the  isolated  point  Pi  on  a  mean  stress  contour  Si 
and  let  a  local  orthonormal  basis  be  fixed  in  P.^  such  that 
the  tangent  and  normal  to  at  P^  are  the  mutually 
perpendicular  axes  u  and  v  as  depicted  in  Fig.  5  (a,b) . 
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A  C.-  VPj. 

AS  *  t'cbnJ. 

SC  -  'Vp  -•«•/»•  fa.nl 

■/  ^  •/ 
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Fig.  5.5a  System  Coordinates  :  XL- positive 
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A3  -  'Vpj-  fan* 

&~f*i  -  AApj.  —  'Vpj  fan* 

DPt'  -  C^tcpj  -  nrPj  tan  *  )  Cos*. 

Qpj  -  $>.•  —(*1  pj  ~ '*Pj  f*n*-)  coS* 


Fig.  5.5b  System  Coordinates  :  AX- negative 


A  neighbouring  contour  Sj  must  touch  the  circle  which  has 
the  tangent  u  at  and  centre  on  the  normal.  The  circle 
will  satisfy  two  of  the  three  constraints  if  it  touches 
but  does  not  intersect  the  curve  S  j .  In  order  to  find  the 
centre  of  this  inscribed  circle  the  length  of  the  radius 
must  be  minimised  over  the  set  of  all  the  circles  that  can 
be  drawn  through  with  a  common  tangent  p  as 
depicted  in  Fig.  5.5a  or  Fig.  5.5b.  Because  the 

set  is  infinite  and  unbounded  it  is  necessary  to  confine 
the  search  to  a  fixed  interval  measured  along  the  curve  Sj. 
Let  the  interval  have  as  end  points  the  intersection  of  the 
normal  from  Pi  on  Sj  and  an  arbitrary  point  denoted  by  Pj 
on  S j .  The  equation  of  the  circles  that  intersect  the 
curve  Sj  in  the  interval  is  given  in  implicit  form  by  the 
expression: - 

u^  -  2vr  +  v^  =  0  5.16 

and  the  centres  are  located  at  u  and  v  of  the  local  co- 

c  c 

ordinate  system  where  it  follows  from  Equation  5.16. 


n  “  1,2,3,... 


5.17 


The  circle  of  radius  r  =  |v  |  is  the  required  member  of 

n  cn 

the  family  when  rn  is  an  extremum  value.  The  minimum  value 
is  easily  found  by  the  Fibonacci  search  method  implemented 
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by  sampling  along  the  curve  in  the  prescribed  interval.* 

The  remaining  constraint  is  satisfied  by  insisting  that 
two  neighbouring  circles  bounded  by  and  have  their 
centres  separated  by  the  sum  of  their  respective  radii 
taken  along  the  curve  joining  the  centres.  To  implement 
this  constraint  let  the  initial  trial  circles  intersect 
i.e.  the  radical  axis  lies  on  both  circles.  In  subsequent 
steps  the  intersection  on  the  classical  Venn-type  diagram 
is  reduced  to  a  null  set  which  then  satisfies  the  third 
constraint. 


Thus  far  the  point  is  arbitrarily  positioned  on  the 
curve  .  In  order  to  locate  the  first  trajectory  to 
we  need  a  symmetry  axis  where  the  circle  that  starts  the 
process  can  be  defined  without  resort  to  the  search 
technique.  The  radius  of  the  initial  circle  in  each  flow 
tube  is  calculated  from  the  intercept  of  S.  and  Sj  on  the 
axis  of  symmetry.  It  transpires  that  a  good  first  approxi¬ 
mation  to  the  location  of  is  given  by  offsetting  PA  a 
distance  of  2r  along  the  curve  S .  where  r  represents 
the  radius  of  the  predefined  circle.  The  point  P^  is 
located  at  a  distance  r  along  S . .  For  either  parallel 
or  diverging  curves  (S.^  and  S..)  this  choice  ensures  that 
the  first  trial  circle  intersects  the  previous  best  fit 
to  the  given  constraints.  The  final  position  of  P^  on 


* 


The  search  employs  the  Fibonacci  sequence  of  integer 
values  given  by: 


Alternatively  a  tabulation  of  sampling  stations  based 
on  this  formula  is  available.  Henley  and  Williams 


(1972) . 
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the  curve  S.^  is  found  by  iteration  using  as  incremental 
steps  the  distances  A  given  by  the  Expression  5.18  and 
repeated  application  of  the  Fibonacci  search, 

.  2rn+l  6  5.18 


where  <5  denotes  the  overlap  of  the  two  circles  alone 
the  curve  through  their  centres  cf.  Fig.  5.6 
rn  and  rn+1  are  the  radii  of  the  neighbouring  circles. 

If  the  curves  S.^  and  derive  from  a  functional  relation¬ 
ship  such  as  the  Eqn.  5.3,  it  merely  remains  to  relate 
the  local  coordinate  system  to  the  global  axes.  For  the 
sake  of  generality  the  point  Pj  is  taken  to  right  of  the 

normal  P  in  addition  to  the  left  of  the  point  which  is  the 
n 

desired  location  of  Pj  in  this  context  as  shown  in  Fig.  5.5.b 
The  coordinate  systems  provide  the  relationships: 

fp  =  -  (up  -  v  tana)  cosa 

3  i  j  j 

fp  -  fp  +  v  seca  +  (u  -  v  tana)sina 

j  *i  3  Fj  j 

5.19 

where  tan  a  is  the  slope  of  the  tangent  at  P^,  and  (4>p,fp) 
represent  global  coordinates  of  points  on  the  mean  stress 
contours.  By  matrix  inversion  the  values  of  up,vp  become 
the  dependent  variables  required  for  evaluating  the  radii 
of  the  inscribed  circles.  Of  course  in  a  computer  program 
it  is  not  necessary  to  graph  the  circles;  the  global  co¬ 
ordinates  of  the  points  of  tangency  and  the  location  of 
centres  is  all  that  is  required  from  the  iteration  process. 
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The  trajectories  representing  the  flow  lines  pass  through 
the  contact  points  of  the  inscribed  circles,  hence  it  is 
necessary  to  generate  a  curve  through  the  centres  of  the 
circles,  and  subsequently  to  find  the  contact  points  by 
offsetting  chord  lengths  along  the  median  curve  equal  to 
the  sum  of  the  radii  of  all  feasible  circles  within  the 
flow  tube. 


At  first  sight  the  procedure  outlined  above  might  seem 
complicated;  in  fact  it  is  a  straightforward  programming 
exercise  with  only  one  exceptional  difficulty.  The 
difficulty  forecasted  is  that  of  finding  the  intersection 
point  of  the  normal  which  bounds  the  search  level.  The 
difficulty  is  exacerbated  by  the  use  of  parametric  cubic 
splines.  The  normal  at  the  point  ( <|>^ ,  4^)  is  given  by 
(cf .  Fig.  5.5) : 

V  -  4^  =  —  <  4>  -  (j^)  cota 

and  the  spline  fit  to  S..  is  of  the  form  (Equation  4.16): 

4  =  a-^t3  +  bjt2  +  cjt  +  dj 

,  =  .  /  +  V2  + 

J 


Thus,  the  coordinates  of  the  intersection  must  be  deduced 
by  solving  for  the  parameter  t  between  those  three 
equations  and  subsequently  finding  the  global  coordinates 
represented  by  t.  Here  the  t-parameter  is  one  of  the 
roots  of  a  cubic  equation  viz. 


3  2 

(aT  cota  +  e_)t  +  (b7  cota  +  f,)t  +  (cT  cota  +  gT)t  + 

J  J  J  JO  0 

k  =  0 


where  k  =  V  .  +  e.cota-  hT  -  dTcota 

XX  U  J 
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5.20 


STEP 


1.  Loading  : 

2.  Stress  Distribution  : 

3.  Equi  Stress  Lines  : 

4.  Interpolation  : 

5.  Trajectories  : 

6.  Mesh  Description  : 

7.  Cell  Dimension  : 

8.  Data  Transfer  : 

5.7  Flow  Chart  :  Numerical  sol 
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OPERATION 


Specify  boundary  conditions 
Approximate  tractions  by 
continuous  functions  on 
boundary. 

Solve  finite  difference 
equations  for  Stress  Components. 

Plot  contours  of  selected 
stress  component  B-lines. 

Pass  spline  function  through 
stress  contours. 

Generate  orthogonal  trajectories 
to  stress  contours  by  search 
routine  a-lines. 

Store  Cartesean  coordinates 
of  a-B  intersections. 

Find  cell  dimensions  and 
Jacobian  transformation  of 
areas . 

Transfer  data  to  flow  Analysis 
Program. 


of  Potential  Stress  Field 


Because  of  multiple  roots  of  which  only  one  is  relevant,  the 
Eqn.  5 . 20  can  prove  a  nuisance  in  numerical  solutions .  Indeed, 
the  main  difficulty  arises  from  the  fact  that  the  constants 
aT,  b,  etc.  are  specific  to  the  knots  on  the  spline  and 
hence  it  is  necessary  to  specify  the  regime  of  parameter 
t  beforehand.  One  can  only  hope  that  this  problem  can  be 
resolved  by  neat  programming!  If  this  drawback  is  removed 
it  may  even  become  practicable  to  find  the  coordinates  of 
the  intersection  of  the  trajectories  which  is  regarded  as 
the  final  (but  not  essential)  step  in  the  proposed  method. 
The  procedure  is  illustrated  by  the  flowchart  in  Fig.  5.7. 
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CHAPTER  6 

Finite  Element  Analyses  of  Creep 
and 

Continuum  Mechanics 


j 


i 


! 


6.1  Introduction 


In  this  Chapter  a  brief  review  of  the  application  of 
finite  elements  is  presented  as  background  to  further 
development  of  the  method  in  this  research  project. 
The  review  is  mainly  prepared  from  the  formulation 
given  by  Penny  and  Marriott  (1971) ;  this  treatment 
is  elementary  but  has  the  appeal  that  it  contains  the 
fundamentals  which  form  the  basis  of  subsequent 
literature  on  the  topic.  The  present  author  has 
endeavoured  to  concentrate  on  those  sections  that  are 
not  held  in  common  by  the  various  contributors  (as 
for  instance  the  peripheral  development  of  solution 
routines  for  time  stepping/instability  caused  by  ill- 
conditioned  matrices  and  like  problems) .  The 
author's  contribution  centres  around  the  manipulation 
of  reference  frames  and  the  handling  of  creep  test 
data  in  analyses. 

6.2  Summary  of  Matrix  Displacement  Method 

The  initial  loading  problem  is  solved  first;  initial 
strains  due  to  plasticity,  creep,  thermal  variations 
are  denoted  by  a  vector  (e  }  and  the  matrix  algebra 

i 

is  in  the  usual  notation.  Superscript  'e'  refers 
to  element  quantities.  Subscript  'el'  refers  to 
elastic  properties  and  subscript  c  to  creep. 

Primed  quantities  refer  to  local  element  coordinates. 

The  procedure  is  as  follows: 
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0 

(a)  Relate  element  strains  {e  }  to  element 

0  I 

nodal  displacements  {6  }  using  a  convenient 

local  set  of  coordinates 

{ee}  =  [bJ  { <5e ' } 

(b)  Relate  element  stresses  {o  }  to  element  strains 

{ae}  =  [D]Uel}  =  fDj  ({ee}-{ea 

(c)  Construct  element  force/displacement  equations 
in  local  coordinates 

{Fe'}  =  [ke ' J { 6e ' }  {FJ} 

where  {F  }  are  the  nodal  forces 

[k&']  =  XolWT  NTs]  d(vol) 

{Fr,}  =  fvoiWT  >  d^vo1) 

(d)  Construct  a  matrix  transforming  local 
coordinates  {x'}  to  global  coordinates  {x}, 
or  alternatively  specify  displacements  in 
terms  of  shape  functions  referred  to  nodal 
coordinates . 

(e)  Solve  assembled  equilibrium  equations  for 

0 

{6  }  and  element  stresses 

{oe}  =  [d]([b]  [A_l}{5e}-{ei }) 
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(f)  Calculate  creep  rates  using  a  creep  law  in 
terms  of  equivalent  stress  components  or 
stress  invariants 


5  ■  il>*-v2+  (w!  + 

+  6(^y  +  Tyz  +  ’zx'f 

(ex-£y) 2  +  (ey"eZ>2  +  (EZ_ex)2 

+  6  ( T  2  +t2  +t2  )"V^ 

xy  yz  zx  1 

(g)  Calculate  the  displacement  and  stress  rates 

with  replacement  of  all  terms  of  the  elastic 

analysis  by  their  rates.  In  particular, 

for  creep  {e  }  becomes  {£  } 

1  c 

Ue]  =  [b] 

{ae}  =  -  U^) 

(h)  Choose  a  small  time  interval  Lx  and  calculate 
quantities  at  the  end  of  that  interval  e.g. 
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{0  W  -  {0  }T  +  {4  }AT 


(i)  Repeat  for  successive  time  increments  until  a 
stationary  state  is  achieved. 

This  summarizes  the  finite  element  technique  if 
a  constitutive  equation  is  adopted  to  describe 
material  behaviour. 


Comment : 


The  formulation  outlined  above  is  particularly 
applicable  to  forced  flow  where  the  stress 
distribution  as  a  function  of  time  is  the 
main  interest. 


For  the  sake  of  completeness  the  equivalent  creep  stress 
vs-strain  relationships  for  a  body  of  revolution  with 
axial  symmetry  are  listed  as  follows: 

* 

Ae  /  a a+a r*\ 

Lzr,c  =  j(°x  '  -V5) 


★ 

Ac  „  o  +a 

T<  6  '  -Vs’ 


At  =  - ( Ae  +e  ) 

z,c  v  r,c  v ,c 


* 

3  A^ 

2  *  Trz 


where 


*  -  ■^({ar-o0)2+(oe“°z)2+(az_ar)2+6TrzI  2 

ie  =  Ae_  “Ae„  _)2+(Ae„  _-Ae_  _)  ^ 


6.3  Discretization  based  on  Stress  Fields 


In  a  continuum  problem  the  field  variable,  whether  it  be 
stress,  displacement  or  some  other  quantity,  possesses 
infinitely  many  values  because  it  is  a  function  of  each 
generic  point  in  the  region.  The  finite  element  process 
reduces  the  problem  to  one  of  a  finite  number  of  unknowns 
at  prescribed  locations  (element  nodes)  in  terms  of  assumed 
approximating  functions  within  each  element.  This  dis¬ 
cretization  known  as  mesh  generation  is  ordinarily  dictated 
by  the  geometrical  configuration  of  the  solution  domain  and 
by  the  implementation  of  the  boundary  conditions.  However, 
the  process  of  discretizing  the  interior  of  the  region  is 
generally  one  of  analytical  preference;  a  mesh  which 
consists  of  triangles  (straight  or  curved  sides)  meets  most 
applications.  It  has  been  the  usual  practice  to  attempt  the 
simultaneous  evaluation  of  stress  and  displacement  in  the 
primary  analysis  using  elaborate  finite  element  programs, 

Key  et  al.  (1978).  This  approach  is  standard  practice  if 
the  stress  field  is  not  invariant  with  the  elapsed  time  of 
straining  the  material.  However,  for  those  problems  where 
the  stresses  in  the  steady  state  have  essentially  the  same 
magnitude  as  the  elastic  stresses  the  development  thus  far 
indicates  that  a  simplistic  approach  is  possible,  which 
obviates  the  need  for  a  finite  element  analysis.  Once  the 
stresses  are  evaluated  it  merely  remains  to  introduce  some 
form  of  constitutive  relationship  to  establish  the  dis¬ 
placement  field  by  integration. 

6.4  Incremental  Solution  of  Displacement  Field 

At  the  outset  the  stress  distribution  on  a  cell  bounded  by 
mean  normal  stress  contours  (isopacics)  and  flow  lines  was 
illustrated  in  Fig.  2. l.For  the  cases  of  plane  strain  and 


( 10  3) 


axisynunetric  problems  the  mean  stress  distributions  are 
determinable  by  the  method  described  in  Chapter  5  (or 
indeed  by  the  more  conventional  method,  Airy  stress 
functions,  boundary  elements  (BIEM)  or  finite  elements 
(FEM)  . 

Assume  that  it  is  possible  to  relate  the  straining  of  the 
material  to  localised  values  of  a  function  of  the  stresses; 
for  instance  the  normal  stress  gradient  is  the  function 
suggested  in  Chapter  2.  On  this  basis  the  displacement 
field  will  be  evaluated  by  using  the  relationships  expressed 
in  the  theory  of  finite  strain  and  the  coordinate  transforma¬ 
tions  of  continuum  mechanics .  Herein  the  pertinent  informa¬ 
tion  follows  closely  the  treatment  of  continuum  mechanics 
to  be  found  in  texts  by  Sokolnikoff  (1967)  and  Hodge  (1970) . 
For  the  sake  of  unification  the  author  has  taken  some 
liberties  with  the  notations  in  these  texts,  and  an  attempt 
has  been  made  to  couch  the  analysis  in  terminology  familiar 
to  readers  of  engineering  literature. 

According  to  Sokolinikoff  the  deformation  of  a  continuum 
medium  can  best  be  described  with  respect  to  three  reference 
frames  as  shown  in  Fig.  6.1;  a  fixed  reference  frame  Y 
determined  by  the  basis  c^  a  moving  reference  frame  c  with 
basis  b^  and  a  fixed  reference  frame  x  with  basis  a^.  The 
undeformed  region  is  denoted  by  and  the  same  region  after 
deformation  at  elapsed  time  t  is  denoted  by 
Two  points  PQ  and  will  arrive  at  the  positions  p  and  p  1 
in  respectively,  if  the  continuum  is  disturbed  by 
external  (or  internal)  tractions.  The  relationship  between 
the  disturbed  and  undisturbed  positions  of  PQ  and  P^  can  be 
established  in  terms  of  the  three  coordinate  systems;  the 
Y  system  is  Cartesian  and  the  X  -  systems  in  general  are 
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curvilinear.*  Reference  to  Fig.  6.1  leads  to  the  following 

?i  =  i  2  3 

(x  ,  /  X3,  t) 

with  inverse 

x1  =  xi(51,C2,?3,  t) 

A  material  point  PQ  in  relative  to  the  orthogonal  Cartesian 
frame  Y  is  determined  by  the  position  vector 

ro  =  Ciyjtx1,  x2,  x3,  t) 


The  location  of  the  same  point  P  in  t  is  determined  by  the 
vector 

r  =  ci  y1*?1,?2,?3,  t) 


The  base  vectors  bj  in  the  moving  frame  £  are  given  by 


K  _  <5r  _  „ 

'  7?  ' Ci 


<sy  (x,t) 
6£J 


The  corresponding  base  vectors  in  are  denoted  by 


a,  =  6ro 
Sx 


C, 

6xJ 


12  3 

The  labels  (x  ,  x  ,  x  )  of  a  given  material  point  P  in  the 

12  3  ° 

X-frame  and  (£  ,  £  ,  £  )  in  the  £-frame  have  the  same  values 

and  transposition  of  the  indices  is  adopted  for  ease  of 

writing  as  the  occasion  presents  (except  for  base  vectors 

for  which  reciprocal  basis  are  a* ,b* , c*) .  With  reference 

to  Fig.  6.1  the  following  relationships  can  be  written  by 

inspection  of  the  initial  and  final  geometry. 


*  For  generality  the  curvilinear  coordinates  are  not 
necessarily  orthogonal  sets. 
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Vector  OP  = 
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l 

or  b^^  = 

a .  + 
i  6x^ 

The  displacement  vector  r  of  the  point  PQ  is  defined  by 

^  .  .Vi 

r  =  r  -  ro 

and  the  component  of  r  relative  to  basis  a^  by  (u,v,w)?  and  its 
components  relative  to  basis  b^  by  (u^v^w1). 

The  vector  ?QPg  =  dr^  which  can  be  expressed  in  the  form 
drQ  =  a^dxi 

and  the  square  of  the  arc  elements  ds  in  n  is 

o  o 

2 

(dso)  =  dVdro  =  aIajdxidxJ 

=  hiJdx.dxJ 

where  h^  denotes  the  metric  coefficient.  6.1 
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Similarly  the  square  of  the  element  of  arc  ds  in  flt  determined 
by  the  corresponding  vector  PP'  can  be  expressed  as 

dr  =  b^dx.^ 

2 

.*.  (ds)  =  b^*bjdx^dXj 

giJ  dxi  dxJ  6*2 

Again  g^  is  the  metric  coefficient  (the  dot  product  of  two 
vectors  and  is  a  scalar) . 

In  expanded  form  this  reads 


(ds)2  = 

gll(dxl)2 

+ 

g12dxldx2 

+ 

g13dxldx3 

+ 

g21dx2dxl 

+ 

g22(dx2)2 

+ 

g23dx2dx3 

+ 

g31dX3dxl 

+ 

g32dx3dx2 

+ 

g33 (dx3^ 

or  in  compact  quadratic  form 


'gn 

g12  g13 

dx1' 

(ds)2 

[dxl 

dx2  dx3l 

g21 

g22  g23 

4 

dx2 

.  g31 

g32  g33 

dx3 

. 

= 

(dXj,) 

(gij)  dx 

J  ' 

wherein 

giJ 

3 

k=l 

iUs 

Sxi 

iJJs 

5xj 

i 

,j  =  1,2,3 
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i.e. 
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2 

2 

55,  2 

gll 

+ 

( — -) 
6xl 

9i2 

Hi 

4- 
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3  <■! 
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Sx.^ 

6x2 

6X^ 
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6 . 5  Strain  Rate 


Finite  strain  is  defined  by  one  of  the  following:- 

(a)  Green's  strain  tensor 

(b)  Almansi's  strain  tensor 

(c)  Natural/logarithmic  strain 


The  Green's  tensor  gives  quadratic  strains  with  respect  to  the 
initial  coordinates  by  the  definition: - 


where  denotes  the  Kronecker  delta  and  ^  is  the  length 

ratio  of  the  line  element  initially  parallel  °  to  the  x^-axis 
which  may  have  any  current  orientation. 


The  A1 man si 


tensor  reads 

6xk  5xk\ 

jtJ 


6.5 


Here  is  the  length  ratio  of  the  element  parallel  to  the 
current  °  axis,  which  generally  had  some  other  orientation  in 
the  initial  state. 


Thus  the  Green's  strain  tensor,  which  is  the  most  commonly  used 
measure  of  finite  strain,  gives 


(ds)2  -  (dsQ)2 


2eiJ  dxi  dxj 


6.6 


=  e the  set  of  functions 


Since  (6.6)  is  an  invariant  and 
e i j ( x , t )  represent  a  tensor  L  with  respect  to  a  class  of 
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admissible  transformation  of  coordinates  X,  with  the  basis 
ai  covering  the  region  fi  .  The  same  set  of  functions 
Gij(C/t)  also  determines  a  tensor  E  with  respect  to  a  set  of 
transformations  of  coordinates  determined  by  the  basis  bi  of 
the  final  state  n.  The  first  interpretation  is  termed 
Lagrangian;  it  leads  to  expressions  for  the  displacement 
from  the  position  in  the  undeformed  state  to  the  deformed 
position  On  the  other  hand  the  second  interpretation, 

the  Eulerian  system,  leads  to  expressions  for  displacements 
that  must  have  taken  place  to  get  to  the  final  position  £ 
from  the  undeformed  configuration  x^.  The  distinction  between 
lagrangian  and  eulerian  variables  becomes  negligible  for  small 
strains  and  rotations. 

The  set  of  differential  equations  for  the  components  of  the 
displacement  when  the  functions  are  specified  becomes 

giJ  "  hiJ  =  2eiJ  2eiJ  6,7 

Therefore,  from  a  purely  geometrical  standpoint,  we  may  bring 
a  neighbourhood  of  P  from  its  initial  to  its  final  state  by 
first  translating  the  neighbourhood  to  move  Pq  to  P,  and  then  by 
rotating  it  so  that  orthogonal  coordinate  axes  of  the  initial 
state  point  in  the  orthogonal  directions  of  the  final  state, 
and  finally  deforming  the  neighbourhood  by  stretching  the 
three  axes  by  the  ratios  of  the  current  to  the  initial  lengths 
parallel  to  the  initial  axes.  Hodge  (1970)  pp.  135-137. 

Hence  an  arbitrary  translation  of  a  deforming  element  will 

derive  from  two  parts  represented  respectively  by  the  strain 

tensor  T  and  a  rotation  tensor  w4T.  The  elements  of  the 
l  J  i  J 

strain  tensor  are,  in  a  convenient  notation,  given  by 


(110) 


iJ  =■ 


C11 

£  12 

c  13  ' 

eiJ 

=  , 
6Xi 

_iV. 

«XJ 

=  J 

21 

22 

23 

t 

_  1/6  u 

+  J_v\ 

'I  ^ 

c  31 

e  32 

G  33 

2[6Xj 

6X±  J 

1  r 

4 

6.8 

tensor  is 

defined  by 
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As  each  line  element  has  its  own  angle  of  rotation  the  entries 
in  w^j  merely  gives  the  average  rotation  components  of  all  the 
line  segments  in  the  elemental  body;  the  components  of  w^ 
represent  rigid  body  rotations,  Dym  and  Shames  (1973),  Boresi/1974). 
Now  let  Pq  (x^,  x2,  x^)  devote  a  point  at  which  displacement 
u°  and  rotation  w°j  are  known  in  a  simply  connected  region. 

Then  the  displacement  at  any  other  point  P  is  obtained  by 
integration  along  a  continuous  curve  from  PQ  to  P.  If  such  a 
line  integral  is  to  be  independent  of  the  path  it  follows  that 
the  strain  distribution  must  result  in  an  integrand  which 
possesses  the  property  of  an  exact  differential  Boc  cit  Indeed 
this  is  the  criterion  used  in  elasticity  theory  to  set  up  the 
compatibility  equations;  it  forms  the  solitary  condition  on 
the  previous  set  of  differential  equations Eqn.  6.7.  Recall 
that  in  Chapter  11  the  issue  of  continuity  arose  in  connection 
with  strain  across  contiguous  flow  paths.  With  the  stipulation 
of  the  exact  differential  the  displacement  of  the  point  P 
becomes  „  _ 


P  P 

ui  =  ui  +  /£iJ  dxJ  +  /w iJ  dxJ 
P 

o  o 


6.10 
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where 


For  the  general  case  of  an  arbitrary  integration  path  the  dis¬ 
placement  field  must  be  determined  by  integrating  the  rather 
cumbersome  integrands  of  Eqnj6-10with  respect  to  non-rectilinear 
coordinates.  Fortunately  the  sets  of  equations  6.7  and  6.10 
assume  quite  simple  forms  when  an  orthogonal  curvilinear 
system  with  origin  on  a  flow  line  is  chosen.*  By  setting  a 
principal  coordinate  direction  tangential  to  the  flow-line  and 
by  making  the  flow  line  a  path  of  integration  the  equation  for 
the  displacements  becomes  (with  reference  to  Fig.  6.1) 

r.  1  t+At 

Up  ( t )  =  u°  +  /**n  (t)  df  J  6.11 

where  ’e^(t)  denotes  the  direct  strain  rate  component  as  a 
function  of  time  of  straining.  On  an  axis  of  symmetry,  or  on 


*  Along  a  flow  line  the  following  relationships  hold: 

1.  rotations  w^  =  0  due  to  orthogonality  of  coordinate  axes 

2.  metric  coefficients  reduce  to  g..,  h„  and  the  quadratic 

form  becomes  11  11 

(ds)2  =  g11(dx1)2  +  g22(dx2)2  +  g3;J(dx3)2 

1-e.  9U  *  921  "  al*a2  =  0  ’31  "  S>13  “  Va3  '  0 

3.  the  displacements  parallel  to  the  flow  lines  are  sufficient 
to  plot  a  flow  field.  Transverse  component  perpendicular 
to  flow  line  is  determinate  by  property  of  exact  differen¬ 
tial  if  required. 
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a  rigid  boundary  u°  is  zero;  the  value  of  up  represents  the 
extension,  or  contraction,  of  a  flow  path  relative  to  its 
length  at  time  t  measured  from  a  fixed  datum,  i.e.  a  symmetry 
axis  or  rigid  boundary. 

Where  the  boundary  itself  is  subject  to  displacements  of 
interest,  as  beneath  a  surface  traction,  the  Equation  6.11 
provides  a  relationship  for  incremental  solution  of  the 
points  of  the  flow  line  with  respect  to  the  initial  position. 
The  displacements  at  points  along  the  flow  line  are  further 
governed  by  the  Equation  2.10  and  2.11.  (There  is  a  direct 
analogy  here  with  the  problem  of  a  non-uniform  variation  of 
temperature  in  a  heated  elastic  solid;  the  strains  at  any 
point  is  given  by  e  =  aT  and  these  must  be  integrated  with 
respect  to  the  coordinates  to  arrive  at  the  displacements) . 

Evidently  the  evaluation  of  the  displacements  along  individual 
flow  lines  will  yield  the  complete  displacement  field  for 
axially  symmetric  and  plane  strain  problems.  The  displacements 
so  determined  over  incremental  time  steps,  assuming  a  linear 
variation  over  the  duration  of  each  time  step,  will  sum  to  the 
total  displacements  at  any  desired  time  t;  i.e.  the  motion 
between  x^  and  x^_^  is  linear.  As  a  consequence,  the  incre¬ 
mental  velocity  given  by  v^+J^  =  Ax1/At  is  constant  over  the 
time  increment.  The  author  predicts  that  the  application  of 
Eqn.  6.11  will  give  an  enhanced  insight  into  flow  mechanisms 
compared  with  solutions  by  the  finite  element  or  kernel 
methods  discussed  earlier  in  this  report. 
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CHAPTER  7 


INTERPRETATION  OF  CREEP  TESTS  AND 
SUMMARY  OF  INVESTIGATION 


7.1  The  Theory  of  Creep  Potential  as  an  aid  to  Interpreta¬ 
tion  of  Experimental  Data. 


The  theory  of  creep  potential  is  based  on  concepts  intro¬ 
duced  in  plasticity  theory  for  the  purpose  of  evaluating 
relative  strain  rates  in  ideal  plastic  solids.  The 
plastic  potential  function  f(a^j)  is  a  scalar  function  of 
stress  which  yields  relative  strain  rates  (the  time  variable 
is  a  dummy  term)  according  to  the  expression  Drucker  and 
Prager  (1952) 


iJ 


=  A 


*f(°iJ> 

5oj 
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where  A  is  a  non  negative  constant.  Extension  of  this 
concept  to  creep  leads  to  an  analogous  relationship.,  Penny 
and  Marriott  (1971). 


de 


iJ,  c 
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For  isotropic  materials  the  experimental  evidence  seems 
to  favour  the  von  Mises  criterion  for  formulation  of  creep 
potential  surfaces.  In  the  case  of  isotropic  multiaxial 
creep  the  creep  potential  deduced  from  the  criterion  is 


Then  the  creep  strain  increments  are  related  to  the  stress 
components  in  the  expression,  loc  cit. 
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where  =  6\p/6o^j  is  termed  the  stress  deviator.  A 
similar  interpretation  is  employed  in  visco-plasticity 
theory  of  Chapter  I  of  this  report.  However,  the  application 
of  creep  potential  theory  has  not  been  extensive  compared 
with  the  hereditary  integral  theory  or  with  the  empirical 
constitutive  relationship  approach. 


The  present  author  suggests  that  the  concept  of  creep 
potential  deserves  more  attention  than  hitherto  found  in  the 
literature.  It  is  one  approach  that  promises  to  inject  some 
unification  into  the  handling  of  experimental  data.  In  order 
to  generalise,  it  is  imperative  that  different  stress  states 
can  be  taken  into  account  and  that  no  matter  what  is  the 
stress  state  in  the  material  the  interpretation  of  experi¬ 
mental  data  leads  to  a  unigue  value  of  the  velocity  at  each 
point  in  the  body.  To  fix  ideas  consider  the  types  of  data 
derived  from  the  following  set  of  laboratory  tests: 

1.  uniaxial  tension  or  compression  on  prismatic  specimens 

2.  triaxial  compression  test  on  cylindrical  specimens 

3.  hollow  cylinder  subjected  to  differential  lateral 
pressures 


In  uniaxial  tension  there  is  no  variation  in  stress  normal  to 
direction  of  flow  hence  the  notion  of  creep  dependence  on 
stress  gradient  is  inapplicable.  Similarly,  in  the  triaxial 
test  the  stress  gradient  between  the  symmetry  axis  and 
boundary  is  not  significant  as  demonstrated  by  finite  element 
analysis;  the  axial  stress  is  virtually  constant  over  the 
length  of  the  specimen  and  radial  and  hoop  stresses  attain 
approximately  0.005  per  cent,  of  the  axial  component. 
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A  state  of  pure  shear  implies  that  normal  stresses  do  not 
exist  on  the  planes  of  maximum  shear  stress.  On  the  other 

hand  the  hollow  cylinder  test  exhibits  creep  under  the 
action  of  stress  gradient  in  the  radial  direction  but  the 
sum  of  the  radial  and  hoop  stresses  is  a  constant  value. 

In  the  light  of  the  earlier  discussion  of  stress  state  in 
the  axisymmetric  flow  of  a  semi-infinite  medium  caused  by 
normal  stresses  on  a  part  of  the  boundary  the  hypothesis 
of  normal  stress  gradient  as  the  primary  source  of  the 
forcing  function  appears  to  present  just  one  further  option 
for  interpreting  test  data.  The  problem  is  how  to  assimi¬ 
late  the  results  of  different  test  configurations  in  a 
unified  presentation.  The  earlier  mentioned  use  of 
effective  stress-strain  is  one  such  attempt  at  unification 
as  indeed  the  choice  of  constants  makes  the  components 
consistent  for  uniaxial  and  multiaxial  tests.* 


Again  the  stress  components  octohedral  normal  (c»oct)  and 
octohedral  shear  (Toct)  play  a  similar  role  as  these 
stress  components  are  related  to  the  stress  invariants 
and  give  a  measure  of  the  hydrostatic  and  the  deviatoric 
components  of  the  state  of  stress  at  a  point.  These 
relations  are 
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*  Effective  stress  in  this  context  is  not  to  be  confused 
with  the  meaning  of  the  word  in  Soil  Mechanics  litera¬ 
ture;  perhaps  the  term  'equivalent  stress'  or  'generalised 
stress'  would  avoid  the  ambiguity. 
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The  spatial  derivative  of  <7^)Ct  is  taken  as  a  measure  of 
the  shearing  intensity  in  a  discretized  zone  according  to 
the  foregoing  hypothesis  of  normal  stress  gradient  as  the 
forcing  function.  The  hypotheses  is  merely  a  restatement 
of  the  case  for  creep  that  the  von  Mises  criterion  in 
plasticity  theory  may  be  extended  to  include  dependence 
on  hydrostatic  stress.  The  yield  surface  will  generate 
a  circular  cone  symmetrical  about  the  hydrostatic  axis. 

The  extended  von  Mises  criterion  takes  the  form,  Harr 
(1966)  pp  169 

(«i  ~  o2)2  +  (o2  ”  °i)2  +  (°i  “  =  k(o1  +  a2  +  O3)2 

7.5 

where  k  is  a  function  of  frictional  resistance  characterised 
by  the  angle  4>  in  granular  media.  The  extended  von  Mises 
criterion  is  treated  in  more  detail  by  Ifedai  (Vol.  1 

pp  226-228) . 

Stress  Manifolds 

From  the  foregoing  discussion  it  follows  that  we  seek  a  set 
of  functions  in  an  n-dimensional  stress  space  which  has  the 
property  of  C* -continuity  (at  least  in  the  domain  of 
interest) .  Mathematically  the  entity  of  the  set  is  termed 
a  real  analytic  manifold;  by  definition  it  consists  of  n 
functions  of  n  variables  which  can  be  expanded  in  a 
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convergent  Taylor  series  in  a  neighbourhood  of  each  point# 
Bruhat  et  al.  (1977)  pp  112.  To  be  of  value  in  the  context 
of  integration  for  displacement  fields  the  tangent  planes 
to  the  manifold  must  give  the  strain  rate  in  the  direction 
of  maximum  velocity;  ideally  the  dip  of  the  tangent  plane 
yields  the  maximum  strain  rate.  The  interpretation  of  test 
data  may  thusly  be  reduced  to  forming  the  Sobolev  space  of 
functions  written  symbolically  as 

Wp  (u)  =  Fn(oij'  aij '  *  n=l,2,3#...  7.6 

Some  applications  of  this  rather  advanced  topic  in 
functional  analysis  is  covered  in  the  Prentice-Hall  Civil 
Engineering  and  Engineering  Mechanics  series#  Oden  1979. 
Suffice  it  to  say  that  the  theory  of  topological  spaces 
provides  a  formalised  approach  to  characterising  a  set  of 
functions  and  to  investigating  distributional  partial 
derivatives  of  the  set. 


7.2  Summary 

The  study  of  creep  undertaken  in  this  report  has  examined 
the  subject  from  a  variety  of  standpoints.  The  first 
chapter  deals  with  the  classical  theories  and  serves  as 
background  to  further  development.  As  with  works  of  an 
original  nature  the  outcome  consists  of  a  mixture  of 
disappointments  and  satisfaction.  On  the  negative  side 
the  postulate  of  a  transform  for  quantifying  strain  rate 
is  defective  in  the  sense  that  the  answer  depends  on  the 
discretization  of  the  solution  domain#  cf.  Chapter  11. 
Another  disappointing  aspect  emerges  in  Chapter  111  wherein 
it  is  shown  that  only  a  limited  set  of  flow  patterns  could 
be  generated  by  the  proposed  harmonic  analysis;  the  set 
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contains  few  patterns  symmetric  about  the  axis  of 
symmetry  of  a  loaded  area  on  a  semi-infinite  half  space. 

Apart  from  these  findings  positive  contributions  may  be 
claimed  as  follows: 

1.  a  framework  is  established  for  predicting  the  dis¬ 
placement  field  as  a  function  of  time  of  straining 
provided  that  strain  rates  can  be  specified  at  sufficient 
points  in  the  solution  domain. 

2.  The  stress  fields  are  determined  by  a  summation 
process  in  a  simple  algorithm. 

3.  The  geometry  of  Laplacian  flow  nets  is  treated  in 
detail  as  are  the  different  types  of  boundary  condition. 

4.  Computer  programs  have  been  prepared  for  implementing 
some  of  the  analytical  techniques. 

5.  The  study  culminates  in  the  suggestion  that  functional 
analysis  offers  an  approach  to  a  unified  interpretation 

of  experimental  data. 

6.  The  main  simplification  results  from  the  proposed  method 
of  direct  integration  along  flow  paths.  This  innovation 
allows  a  free  choice  of  the  strain  rate  tensor  because 
infinitesimal  and  finite  strains  can  be  handled  with  equal 
ease  in  evaluating  the  line  integral  for  displacement; 

the  strain  rates  are  formulated  as  a  part  of  the  inter¬ 
pretation  of  experimental  data  with  no  further  implications 
in  the  analytical  technique. 
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